From 7cfdb16e052b10c495eb9143b46d33ac932e4cec Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 23 Sep 2024 12:00:20 -0700 Subject: [PATCH 01/14] Add adjoint load for the case of adjoints associated with resultant/dual forces. --- src/serac/physics/base_physics.hpp | 8 +-- src/serac/physics/solid_mechanics.hpp | 61 ++++++++++--------- .../physics/tests/solid_reaction_adjoint.cpp | 5 +- 3 files changed, 38 insertions(+), 36 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 20690107a..e1e69720d 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -345,10 +345,10 @@ class BasePhysics { * @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) + virtual void assembleDualAdjointLoad(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", + SLIC_ERROR_ROOT(axom::fmt::format("assembleDualAdjointLoad not enabled in physics module {}, dual name {} requested", name_, dual_name)); } @@ -361,7 +361,7 @@ class BasePhysics { * @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 + * @pre `assembleDualAdjointLoad' for the desired dual (reaction) and `reverseAdjointTimestep` must be called before * this */ virtual const serac::FiniteElementDual& computeDualSensitivity(const serac::FiniteElementState& reaction_direction, @@ -379,7 +379,7 @@ class BasePhysics { * qoi * @return reaction sensitivity field * - * @pre `computeDualAdjointLoad' for the desired dual (reaction) and `reverseAdjointTimestep` must be called before + * @pre `assembleDualAdjointLoad' for the desired dual (reaction) and `reverseAdjointTimestep` must be called before * this */ virtual const serac::FiniteElementDual& computeDualShapeSensitivity( diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 3ba8818ab..2080fa922 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1337,6 +1337,33 @@ class SolidMechanics, std::integer_se } } + /// @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_; + }; + /// @overload void reverseAdjointTimestep() override { @@ -1463,7 +1490,7 @@ class SolidMechanics, std::integer_se const serac::FiniteElementDual& reactions() const { return reactions_; }; /// @overload - void computeDualAdjointLoad(const std::string& dual_name, + void assembleDualAdjointLoad(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"); @@ -1471,38 +1498,12 @@ class SolidMechanics, std::integer_se 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; + + // this is currently hard coded for the case of quasi-statics jacobian->MultTranspose(reaction_direction, reactions_adjoint_load_); - setAdjointLoad({{"displacement", reactions_adjoint_load_}}); + displacement_adjoint_load_ -= 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; diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index 2c03485d9..00a2a7db5 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -121,7 +121,9 @@ auto computeSolidMechanicsQoiSensitivities(BasePhysics& solid_solver) auto reactionDirections = createReactionDirection(solid_solver, 0); - solid_solver.computeDualAdjointLoad(solid_solver.dualNames()[0], reactionDirections); + FiniteElementDual zero(solid_solver.state("displacement").space(), "zero"); + solid_solver.setAdjointLoad({{"displacement", zero}}); + solid_solver.assembleDualAdjointLoad(solid_solver.dualNames()[0], reactionDirections); solid_solver.reverseAdjointTimestep(); shear_modulus_sensitivity += solid_solver.computeTimestepSensitivity(0); @@ -177,7 +179,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); } From a5d237b1b64f42f9911e265408a304bc438de25e Mon Sep 17 00:00:00 2001 From: Agent Style Date: Mon, 23 Sep 2024 15:27:25 -0700 Subject: [PATCH 02/14] Apply style updates --- src/serac/physics/base_physics.hpp | 7 ++++--- src/serac/physics/solid_mechanics.hpp | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index e1e69720d..dd2ea506d 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -345,11 +345,12 @@ class BasePhysics { * @param reaction_direction A FiniteElementState which specifies how the reactions dofs are weighted for the reaction * qoi */ - virtual void assembleDualAdjointLoad(const std::string& dual_name, const serac::FiniteElementState& reaction_direction) + virtual void assembleDualAdjointLoad(const std::string& dual_name, + const serac::FiniteElementState& reaction_direction) { (void)reaction_direction; - SLIC_ERROR_ROOT(axom::fmt::format("assembleDualAdjointLoad not enabled in physics module {}, dual name {} requested", - name_, dual_name)); + SLIC_ERROR_ROOT(axom::fmt::format( + "assembleDualAdjointLoad not enabled in physics module {}, dual name {} requested", name_, dual_name)); } /** diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 2080fa922..3add656d8 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1491,7 +1491,7 @@ class SolidMechanics, std::integer_se /// @overload void assembleDualAdjointLoad(const std::string& dual_name, - const serac::FiniteElementState& reaction_direction) override + const serac::FiniteElementState& reaction_direction) override { SLIC_ERROR_ROOT_IF(dual_name != "reactions", "Solid mechanics has reactions as its only dual"); From 296166c5ec25decf7aeb68a8b2eddeb8fc3b15d3 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 3 Oct 2024 12:19:46 -0700 Subject: [PATCH 03/14] Add test that combined reaction and displacement adjoint gives correct answers. --- .../tests/quasistatic_solid_adjoint.cpp | 8 ++--- .../physics/tests/solid_reaction_adjoint.cpp | 35 ++++++++++--------- 2 files changed, 22 insertions(+), 21 deletions(-) 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 00a2a7db5..d0cee923a 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 { @@ -29,8 +30,9 @@ using SolidMechanicsType = SolidMechanics, H1<1>>>; 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::ParameterizedNeoHookeanSolid; +using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; +auto geoNonlinear = GeometricNonlinearities::Off; constexpr double boundary_disp = 0.013; constexpr double shear_modulus_value = 1.0; @@ -41,7 +43,7 @@ std::unique_ptr createNonlinearSolidMechanicsSolver(mfem::Pa const SolidMaterial& mat) { static int iter = 0; - auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, + auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, solid_mechanics::default_quasistatic_options, geoNonlinear, physics_prefix + std::to_string(iter++), mesh_tag, std::vector{"shear modulus", "bulk modulus"}); @@ -101,12 +103,14 @@ 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"); - 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) @@ -121,8 +125,11 @@ auto computeSolidMechanicsQoiSensitivities(BasePhysics& solid_solver) auto reactionDirections = createReactionDirection(solid_solver, 0); - FiniteElementDual zero(solid_solver.state("displacement").space(), "zero"); - solid_solver.setAdjointLoad({{"displacement", zero}}); + 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.setAdjointLoad({{"displacement", displacement_adjoint_load}}); solid_solver.assembleDualAdjointLoad(solid_solver.dualNames()[0], reactionDirections); solid_solver.reverseAdjointTimestep(); @@ -208,7 +215,7 @@ 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(qoi_base, -0.35, 1e-14); EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); } @@ -224,7 +231,7 @@ 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(qoi_base, -0.35, 1e-14); EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); } @@ -233,12 +240,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; } From e0cf1f83ec2f0336ebe5c6d60559f3491f629128 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 3 Oct 2024 12:20:16 -0700 Subject: [PATCH 04/14] Remove exactnes test which is no longer appropriate. --- src/serac/physics/tests/solid_reaction_adjoint.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index d0cee923a..4548e2ec9 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -215,7 +215,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); } @@ -231,7 +230,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); } From 2126672b9443d795ace64b78e4a6d05d73a0476f Mon Sep 17 00:00:00 2001 From: Agent Style Date: Thu, 3 Oct 2024 12:22:48 -0700 Subject: [PATCH 05/14] Apply style updates --- src/serac/physics/tests/solid_reaction_adjoint.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index 4548e2ec9..fc077b4c7 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -30,7 +30,7 @@ using SolidMechanicsType = SolidMechanics, H1<1>>>; const std::string mesh_tag = "mesh"; const std::string physics_prefix = "solid"; -//using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; +// using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; auto geoNonlinear = GeometricNonlinearities::Off; @@ -43,7 +43,7 @@ std::unique_ptr createNonlinearSolidMechanicsSolver(mfem::Pa const SolidMaterial& mat) { static int iter = 0; - auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, + auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, solid_mechanics::default_quasistatic_options, geoNonlinear, physics_prefix + std::to_string(iter++), mesh_tag, std::vector{"shear modulus", "bulk modulus"}); @@ -106,8 +106,8 @@ double computeSolidMechanicsQoi(BasePhysics& solid_solver) solid_solver.resetStates(); solid_solver.advanceTimestep(0.1); - const FiniteElementDual& reactions = solid_solver.dual("reactions"); - auto reactionDirections = createReactionDirection(solid_solver, 0); + const FiniteElementDual& reactions = solid_solver.dual("reactions"); + auto reactionDirections = createReactionDirection(solid_solver, 0); const FiniteElementState& displacements = solid_solver.state("displacement"); return innerProduct(reactions, reactionDirections) + 0.05 * innerProduct(displacements, displacements); From 06d2f03e8d2bc674610a57737754c596eca415bc Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 7 Oct 2024 16:27:32 -0700 Subject: [PATCH 06/14] Significantly simplify dual adjoint calculations. --- src/serac/physics/base_physics.hpp | 65 +++-------------- src/serac/physics/solid_mechanics.hpp | 70 ++++++++----------- .../physics/tests/solid_reaction_adjoint.cpp | 12 ++-- 3 files changed, 46 insertions(+), 101 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index dd2ea506d..288cbe619 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 assembleDualAdjointLoad(const std::string& dual_name, - const serac::FiniteElementState& reaction_direction) - { - (void)reaction_direction; - SLIC_ERROR_ROOT(axom::fmt::format( - "assembleDualAdjointLoad 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 `assembleDualAdjointLoad' 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 `assembleDualAdjointLoad' 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 @@ -422,6 +370,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 setDualAdjointLoad(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 @@ -551,7 +508,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 3add656d8..9aee24ef3 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -386,6 +386,7 @@ class SolidMechanics, std::integer_se implicit_sensitivity_velocity_start_of_step_ = 0.0; reactions_ = 0.0; + reactions_adjoint_load_ = 0.0; u_ = 0.0; v_ = 0.0; @@ -753,7 +754,7 @@ 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_; @@ -1337,32 +1338,29 @@ class SolidMechanics, std::integer_se } } - /// @overload - const serac::FiniteElementDual& computeDualSensitivity(const serac::FiniteElementState& reaction_direction, - size_t parameter_field) override + virtual void setDualAdjointLoad(std::unordered_map loads) override { - SLIC_ASSERT_MSG(parameter_field < sizeof...(parameter_indices), - axom::fmt::format("Invalid parameter index '{}' requested for reaction sensitivity.")); + SLIC_ERROR_ROOT_IF(loads.size() == 0, "Adjoint load container size must be greater than 0 in the solid mechanics."); - auto drdparam = serac::get(d_residual_d_[parameter_field](time_end_step_)); - auto drdparam_mat = assemble(drdparam); + auto reaction_adjoint_load = loads.find("reactions"); - drdparam_mat->MultTranspose(reaction_direction, *parameters_[parameter_field].sensitivity); + SLIC_ERROR_ROOT_IF(reaction_adjoint_load == loads.end(), "Adjoint load for \"reaction\" not found."); - return *parameters_[parameter_field].sensitivity; - }; + reactions_adjoint_load_ = reaction_adjoint_load->second; + // Add the sign correction to move the term to the RHS + reactions_adjoint_load_ *= 1.0; - /// @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_; - }; + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, + *parameters_[parameter_indices].state...); + std::unique_ptr jacobian = assemble(drdu); + + serac::FiniteElementState displacement_adjoint_load_contribution(displacement_adjoint_load_.space(), "adjoint_load_contribution" ); + + // this is currently hard coded for the case of quasi-statics + jacobian->MultTranspose(reactions_adjoint_load_, displacement_adjoint_load_contribution); + displacement_adjoint_load_ -= displacement_adjoint_load_contribution; + + } /// @overload void reverseAdjointTimestep() override @@ -1439,7 +1437,10 @@ class SolidMechanics, std::integer_se auto drdparam = serac::get(d_residual_d_[parameter_field](time_end_step_)); auto drdparam_mat = assemble(drdparam); - drdparam_mat->MultTranspose(adjoint_displacement_, *parameters_[parameter_field].sensitivity); + serac::FiniteElementState total_adjoint_displacement(adjoint_displacement_); + total_adjoint_displacement += reactions_adjoint_load_; + + drdparam_mat->MultTranspose(total_adjoint_displacement, *parameters_[parameter_field].sensitivity); return *parameters_[parameter_field].sensitivity; } @@ -1453,7 +1454,10 @@ class SolidMechanics, std::integer_se auto drdshape_mat = assemble(drdshape); - drdshape_mat->MultTranspose(adjoint_displacement_, *shape_displacement_sensitivity_); + serac::FiniteElementState total_adjoint_displacement(adjoint_displacement_); + total_adjoint_displacement += reactions_adjoint_load_; + + drdshape_mat->MultTranspose(total_adjoint_displacement, *shape_displacement_sensitivity_); return *shape_displacement_sensitivity_; } @@ -1489,21 +1493,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 assembleDualAdjointLoad(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); - - // this is currently hard coded for the case of quasi-statics - jacobian->MultTranspose(reaction_direction, reactions_adjoint_load_); - displacement_adjoint_load_ -= reactions_adjoint_load_; - } - protected: /// The compile-time finite element trial space for displacement and velocity (H1 of order p) using trial = H1; @@ -1528,6 +1517,7 @@ class SolidMechanics, std::integer_se /// The displacement finite element adjoint state FiniteElementState adjoint_displacement_; + /// The adjoint load (RHS) for the displacement adjoint system solve (downstream -dQOI/d displacement) FiniteElementDual displacement_adjoint_load_; @@ -1547,7 +1537,7 @@ class SolidMechanics, std::integer_se FiniteElementDual reactions_; /// sensitivity of qoi with respect to reaction forces - FiniteElementDual reactions_adjoint_load_; + FiniteElementState reactions_adjoint_load_; /// serac::Functional that is used to calculate the residual and its derivatives std::unique_ptr> residual_; diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index fc077b4c7..1c6be6c9e 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -30,8 +30,8 @@ using SolidMechanicsType = SolidMechanics, H1<1>>>; const std::string mesh_tag = "mesh"; const std::string physics_prefix = "solid"; -// using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; -using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; +using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; +// using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; auto geoNonlinear = GeometricNonlinearities::Off; constexpr double boundary_disp = 0.013; @@ -108,6 +108,7 @@ double computeSolidMechanicsQoi(BasePhysics& solid_solver) const FiniteElementDual& reactions = solid_solver.dual("reactions"); auto reactionDirections = createReactionDirection(solid_solver, 0); + //reactionDirections = solid_solver.dual("reactions"); const FiniteElementState& displacements = solid_solver.state("displacement"); return innerProduct(reactions, reactionDirections) + 0.05 * innerProduct(displacements, displacements); @@ -123,21 +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.setAdjointLoad({{"displacement", displacement_adjoint_load}}); - solid_solver.assembleDualAdjointLoad(solid_solver.dualNames()[0], reactionDirections); + solid_solver.setDualAdjointLoad({{"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); } From 55c2b8a6750cade3b9bf611b0018065ae7a6f818 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 07:15:04 -0700 Subject: [PATCH 07/14] Simplify reaction adjoint implementation even more. --- src/serac/physics/solid_mechanics.hpp | 29 ++++++++------------------- 1 file changed, 8 insertions(+), 21 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 9aee24ef3..209b00fcb 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1347,19 +1347,6 @@ class SolidMechanics, std::integer_se SLIC_ERROR_ROOT_IF(reaction_adjoint_load == loads.end(), "Adjoint load for \"reaction\" not found."); reactions_adjoint_load_ = reaction_adjoint_load->second; - // Add the sign correction to move the term to the RHS - reactions_adjoint_load_ *= 1.0; - - auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].state...); - std::unique_ptr jacobian = assemble(drdu); - - serac::FiniteElementState displacement_adjoint_load_contribution(displacement_adjoint_load_.space(), "adjoint_load_contribution" ); - - // this is currently hard coded for the case of quasi-statics - jacobian->MultTranspose(reactions_adjoint_load_, displacement_adjoint_load_contribution); - displacement_adjoint_load_ -= displacement_adjoint_load_contribution; - } /// @overload @@ -1389,6 +1376,10 @@ class SolidMechanics, std::integer_se auto jacobian = assemble(drdu); auto J_T = std::unique_ptr(jacobian->Transpose()); + serac::FiniteElementState displacement_adjoint_load_contribution(displacement_adjoint_load_.space(), "adjoint_load_contribution" ); + J_T->Mult(reactions_adjoint_load_, displacement_adjoint_load_contribution); + displacement_adjoint_load_ -= displacement_adjoint_load_contribution; + for (const auto& bc : bcs_.essentials()) { bc.apply(*J_T, displacement_adjoint_load_, adjoint_essential); } @@ -1396,6 +1387,8 @@ class SolidMechanics, std::integer_se lin_solver.SetOperator(*J_T); lin_solver.Mult(displacement_adjoint_load_, adjoint_displacement_); + adjoint_displacement_ += reactions_adjoint_load_; + // Reset the equation solver to use the full nonlinear residual operator. MRT, is this needed? nonlin_solver_->setOperator(*residual_with_bcs_); } else { @@ -1437,10 +1430,7 @@ class SolidMechanics, std::integer_se auto drdparam = serac::get(d_residual_d_[parameter_field](time_end_step_)); auto drdparam_mat = assemble(drdparam); - serac::FiniteElementState total_adjoint_displacement(adjoint_displacement_); - total_adjoint_displacement += reactions_adjoint_load_; - - drdparam_mat->MultTranspose(total_adjoint_displacement, *parameters_[parameter_field].sensitivity); + drdparam_mat->MultTranspose(adjoint_displacement_, *parameters_[parameter_field].sensitivity); return *parameters_[parameter_field].sensitivity; } @@ -1454,10 +1444,7 @@ class SolidMechanics, std::integer_se auto drdshape_mat = assemble(drdshape); - serac::FiniteElementState total_adjoint_displacement(adjoint_displacement_); - total_adjoint_displacement += reactions_adjoint_load_; - - drdshape_mat->MultTranspose(total_adjoint_displacement, *shape_displacement_sensitivity_); + drdshape_mat->MultTranspose(adjoint_displacement_, *shape_displacement_sensitivity_); return *shape_displacement_sensitivity_; } From ab4ea4101913555a3531ade8daa8890d988714db Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 07:25:48 -0700 Subject: [PATCH 08/14] Even more simplification of reaction adjoint loads, rename to setDualAdjointBcs to be a bit more mathly accurate. --- src/serac/physics/base_physics.hpp | 2 +- src/serac/physics/solid_mechanics.hpp | 33 +++++++++---------- .../physics/tests/solid_reaction_adjoint.cpp | 2 +- 3 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 288cbe619..82299e151 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -374,7 +374,7 @@ class BasePhysics { * @brief Set the dual loads (dirichlet values) for the adjoint reverse timestep solve * This must be called after setAdjointLoad */ - virtual void setDualAdjointLoad(std::unordered_map) + virtual void setDualAdjointBcs(std::unordered_map) { SLIC_ERROR_ROOT(axom::fmt::format("Adjoint analysis not defined for physics module {}", name_)); } diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 209b00fcb..aa8961b82 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(); @@ -386,7 +386,7 @@ class SolidMechanics, std::integer_se implicit_sensitivity_velocity_start_of_step_ = 0.0; reactions_ = 0.0; - reactions_adjoint_load_ = 0.0; + reactions_adjoint_bcs_ = 0.0; u_ = 0.0; v_ = 0.0; @@ -757,12 +757,12 @@ class SolidMechanics, std::integer_se 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_; } /** @@ -1338,7 +1338,7 @@ class SolidMechanics, std::integer_se } } - virtual void setDualAdjointLoad(std::unordered_map loads) override + virtual void setDualAdjointBcs(std::unordered_map loads) override { SLIC_ERROR_ROOT_IF(loads.size() == 0, "Adjoint load container size must be greater than 0 in the solid mechanics."); @@ -1346,7 +1346,7 @@ class SolidMechanics, std::integer_se SLIC_ERROR_ROOT_IF(reaction_adjoint_load == loads.end(), "Adjoint load for \"reaction\" not found."); - reactions_adjoint_load_ = reaction_adjoint_load->second; + reactions_adjoint_bcs_ = reaction_adjoint_load->second; } /// @overload @@ -1373,22 +1373,21 @@ 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_ = assemble(drdu); + auto J_T = std::unique_ptr(J_->Transpose()); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T); - serac::FiniteElementState displacement_adjoint_load_contribution(displacement_adjoint_load_.space(), "adjoint_load_contribution" ); - J_T->Mult(reactions_adjoint_load_, displacement_adjoint_load_contribution); - displacement_adjoint_load_ -= displacement_adjoint_load_contribution; + 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); lin_solver.Mult(displacement_adjoint_load_, adjoint_displacement_); - adjoint_displacement_ += reactions_adjoint_load_; - // Reset the equation solver to use the full nonlinear residual operator. MRT, is this needed? nonlin_solver_->setOperator(*residual_with_bcs_); } else { @@ -1524,7 +1523,7 @@ class SolidMechanics, std::integer_se FiniteElementDual reactions_; /// sensitivity of qoi with respect to reaction forces - FiniteElementState reactions_adjoint_load_; + FiniteElementState reactions_adjoint_bcs_; /// serac::Functional that is used to calculate the residual and its derivatives std::unique_ptr> residual_; diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index 1c6be6c9e..3c720854e 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -131,7 +131,7 @@ auto computeSolidMechanicsQoiSensitivities(BasePhysics& solid_solver) displacement_adjoint_load *= 0.1; solid_solver.setAdjointLoad({{"displacement", displacement_adjoint_load}}); - solid_solver.setDualAdjointLoad({{"reactions", reaction_adjoint_load}}); + solid_solver.setDualAdjointBcs({{"reactions", reaction_adjoint_load}}); solid_solver.reverseAdjointTimestep(); shear_modulus_sensitivity += solid_solver.computeTimestepSensitivity(0); From d97afa9eae10479c0e9e3153eb9e9579400460e4 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 10:47:14 -0700 Subject: [PATCH 09/14] Addition cleanup, reset sparse matrices to avoid highwater memory. --- src/serac/physics/solid_mechanics.hpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index aa8961b82..302d59b30 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1135,7 +1135,9 @@ class SolidMechanics, std::integer_se SERAC_MARK_FUNCTION; auto [r, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(u), acceleration_, *parameters_[parameter_indices].state...); + J_.reset(); J_ = assemble(drdu); + J_e_.reset(); J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); return *J_; }); @@ -1202,7 +1204,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_; @@ -1354,10 +1358,6 @@ class SolidMechanics, std::integer_se { 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"); @@ -1373,9 +1373,13 @@ class SolidMechanics, std::integer_se if (is_quasistatic_) { auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].state...); - J_ = assemble(drdu); - auto J_T = std::unique_ptr(J_->Transpose()); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T); + 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(); @@ -1413,7 +1417,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_; @@ -1725,7 +1729,9 @@ 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_.reset(); J_ = assemble(drdu); + J_e_.reset(); J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); r *= -1.0; From a43765cbad6a322eef7d3439af4c9da27199a528 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 14:10:37 -0700 Subject: [PATCH 10/14] Test that reaction forces in qois are working when nonlinear and transient. --- src/serac/physics/base_physics.hpp | 13 +++++++ src/serac/physics/solid_mechanics.hpp | 36 ++++++++++++++++--- .../physics/tests/dynamic_solid_adjoint.cpp | 35 +++++++++++------- 3 files changed, 67 insertions(+), 17 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 82299e151..d76973099 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -407,6 +407,19 @@ class BasePhysics { */ 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 cycle The cycle to retrieve state from + * @param state_name The name of the state to retrieve (e.g. "reaction") + * @return The named Finite Element Dual + */ + virtual FiniteElementDual loadCheckpointedDual(const std::string&, int) 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 diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 302d59b30..9322160b6 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -765,6 +765,30 @@ class SolidMechanics, std::integer_se 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_; + } + /** * @brief register a custom domain integral calculation as part of the residual * @@ -1321,9 +1345,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"); @@ -1350,7 +1376,9 @@ class SolidMechanics, std::integer_se SLIC_ERROR_ROOT_IF(reaction_adjoint_load == loads.end(), "Adjoint load for \"reaction\" not found."); - reactions_adjoint_bcs_ = reaction_adjoint_load->second; + if (reaction_adjoint_load != loads.end()) { + reactions_adjoint_bcs_ = reaction_adjoint_load->second; + } } /// @overload diff --git a/src/serac/physics/tests/dynamic_solid_adjoint.cpp b/src/serac/physics/tests/dynamic_solid_adjoint.cpp index fdb2ecda1..075062c55 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,23 @@ 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 +114,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"); + FiniteElementDual reactionsForObjective = solid_solver.dual("reactions"); + double qoi = computeStepQoi(dispForObjective, reactionsForObjective, 0.5 * (dts(0) + dts(1))); - double qoi = computeStepQoi(dispForObjective, 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))); + reactionsForObjective = solid_solver.dual("reactions"); + qoi += computeStepQoi(dispForObjective, reactionsForObjective, 0.5 * (dts(i) + dts(i + 1))); } return qoi; } @@ -138,17 +146,22 @@ std::tuple comp shape_sensitivity = 0.0; 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; } From 396d791a082d535b018a5407be65fb264ff76315 Mon Sep 17 00:00:00 2001 From: Agent Style Date: Tue, 8 Oct 2024 14:12:41 -0700 Subject: [PATCH 11/14] Apply style updates --- src/serac/physics/base_physics.hpp | 4 ++-- src/serac/physics/solid_mechanics.hpp | 24 +++++++++---------- .../physics/tests/dynamic_solid_adjoint.cpp | 22 ++++++++--------- .../physics/tests/solid_reaction_adjoint.cpp | 4 ++-- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index d76973099..e9d936632 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -338,7 +338,6 @@ class BasePhysics { return *shape_displacement_sensitivity_; } - /** * @brief Compute the implicit sensitivity of the quantity of interest with respect to the initial condition fields * @@ -415,7 +414,8 @@ class BasePhysics { * @param state_name The name of the state to retrieve (e.g. "reaction") * @return The named Finite Element Dual */ - virtual FiniteElementDual loadCheckpointedDual(const std::string&, int) const{ + virtual FiniteElementDual loadCheckpointedDual(const std::string&, int) const + { SLIC_ERROR_ROOT(axom::fmt::format("loadCheckpointedDual not enabled in physics module {}", name_)); return *duals_[0]; } diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 9322160b6..4f3a35073 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -385,7 +385,7 @@ 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; @@ -777,15 +777,16 @@ class SolidMechanics, std::integer_se 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...); + 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_)); + SLIC_ERROR_ROOT( + axom::fmt::format("loadCheckpointedDual '{}' requested from solid mechanics module '{}', but it doesn't exist", + dual_name, name_)); return reactions_; } @@ -1160,9 +1161,9 @@ class SolidMechanics, std::integer_se auto [r, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(u), acceleration_, *parameters_[parameter_indices].state...); J_.reset(); - J_ = assemble(drdu); + J_ = assemble(drdu); J_e_.reset(); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); return *J_; }); } @@ -1413,7 +1414,7 @@ class SolidMechanics, std::integer_se 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]; + int j = constrained_dofs[i]; displacement_adjoint_load_[j] = reactions_adjoint_bcs_[j]; } @@ -1535,7 +1536,6 @@ class SolidMechanics, std::integer_se /// The displacement finite element adjoint state FiniteElementState adjoint_displacement_; - /// The adjoint load (RHS) for the displacement adjoint system solve (downstream -dQOI/d displacement) FiniteElementDual displacement_adjoint_load_; @@ -1758,9 +1758,9 @@ class SolidMechanics, std::integer_se auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].previous_state...); J_.reset(); - J_ = assemble(drdu); + J_ = assemble(drdu); J_e_.reset(); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); + 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 075062c55..910ca3380 100644 --- a/src/serac/physics/tests/dynamic_solid_adjoint.cpp +++ b/src/serac/physics/tests/dynamic_solid_adjoint.cpp @@ -49,10 +49,11 @@ double computeStepQoi(const FiniteElementState& displacement, const FiniteElemen FiniteElementState displacement_error(displacement); displacement_error = disp_target; displacement_error.Add(1.0, displacement); - return 0.5 * dt * innerProduct(displacement_error, displacement_error) + 0.05 * dt * innerProduct(reactions, reactions); + return 0.5 * dt * innerProduct(displacement_error, displacement_error) + + 0.05 * dt * innerProduct(reactions, reactions); } -void computeStepAdjointLoad(const FiniteElementState& displacement, const FiniteElementDual& reactions, +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; @@ -114,16 +115,16 @@ 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"); - FiniteElementDual reactionsForObjective = solid_solver.dual("reactions"); - double qoi = computeStepQoi(dispForObjective, reactionsForObjective, 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"); + + dispForObjective = solid_solver.state("displacement"); reactionsForObjective = solid_solver.dual("reactions"); qoi += computeStepQoi(dispForObjective, reactionsForObjective, 0.5 * (dts(i) + dts(i + 1))); } @@ -145,14 +146,14 @@ 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()); + auto previous_reactions = solid_solver.loadCheckpointedDual("reactions", solid_solver.cycle()); computeStepAdjointLoad( previous_displacement, previous_reactions, adjoint_load, adjoint_bcs, 0.5 * (solid_solver.getCheckpointedTimestep(i - 1) + solid_solver.getCheckpointedTimestep(i))); @@ -160,7 +161,6 @@ std::tuple comp 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()); diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index 3c720854e..b154be211 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -32,7 +32,7 @@ const std::string physics_prefix = "solid"; using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; // using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; -auto geoNonlinear = GeometricNonlinearities::Off; +auto geoNonlinear = GeometricNonlinearities::Off; constexpr double boundary_disp = 0.013; constexpr double shear_modulus_value = 1.0; @@ -108,7 +108,7 @@ double computeSolidMechanicsQoi(BasePhysics& solid_solver) const FiniteElementDual& reactions = solid_solver.dual("reactions"); auto reactionDirections = createReactionDirection(solid_solver, 0); - //reactionDirections = solid_solver.dual("reactions"); + // reactionDirections = solid_solver.dual("reactions"); const FiniteElementState& displacements = solid_solver.state("displacement"); return innerProduct(reactions, reactionDirections) + 0.05 * innerProduct(displacements, displacements); From 8efdcecdcf9062cf33dd43f1b4a56f2411f59dae Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 14:20:18 -0700 Subject: [PATCH 12/14] Small naming fix. --- src/serac/physics/solid_mechanics.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 4f3a35073..acb9c281f 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1369,15 +1369,15 @@ class SolidMechanics, std::integer_se } } - virtual void setDualAdjointBcs(std::unordered_map loads) override + virtual void setDualAdjointBcs(std::unordered_map bcs) override { - SLIC_ERROR_ROOT_IF(loads.size() == 0, "Adjoint load container size must be greater than 0 in the solid mechanics."); + SLIC_ERROR_ROOT_IF(bcs.size() == 0, "Adjoint load container size must be greater than 0 in the solid mechanics."); - auto reaction_adjoint_load = loads.find("reactions"); + auto reaction_adjoint_load = bcs.find("reactions"); - SLIC_ERROR_ROOT_IF(reaction_adjoint_load == loads.end(), "Adjoint load for \"reaction\" not found."); + SLIC_ERROR_ROOT_IF(reaction_adjoint_load == bcs.end(), "Adjoint load for \"reaction\" not found."); - if (reaction_adjoint_load != loads.end()) { + if (reaction_adjoint_load != bcs.end()) { reactions_adjoint_bcs_ = reaction_adjoint_load->second; } } From 6ed26e954b6b7231fddc74002fd3d93fd18fbe02 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 8 Oct 2024 16:06:10 -0700 Subject: [PATCH 13/14] Fix docs. --- src/serac/physics/base_physics.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index e9d936632..2d2d54cac 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -400,8 +400,8 @@ 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; @@ -410,11 +410,11 @@ class BasePhysics { * @brief Accessor for getting a single named finite element dual 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. "reaction") + * @param cycle The cycle to retrieve state from * @return The named Finite Element Dual */ - virtual FiniteElementDual loadCheckpointedDual(const std::string&, int) const + 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]; From 1a671414afe7055147a99636560071a70f0cd55f Mon Sep 17 00:00:00 2001 From: Agent Style Date: Tue, 8 Oct 2024 16:15:14 -0700 Subject: [PATCH 14/14] Apply style updates --- src/serac/physics/base_physics.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 2d2d54cac..3d15c57ec 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -414,7 +414,8 @@ class BasePhysics { * @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 + 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];