Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add adjoint load for the case of adjoints associated with resultant/dual forces. #1238

Merged
merged 17 commits into from
Oct 25, 2024
Merged

Add adjoint load for the case of adjoints associated with resultant/dual forces. #1238

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
7cfdb16
Add adjoint load for the case of adjoints associated with resultant/d…
tupek2 Sep 23, 2024
a5d237b
Apply style updates
Sep 23, 2024
296166c
Add test that combined reaction and displacement adjoint gives correc…
tupek2 Oct 3, 2024
e0cf1f8
Remove exactnes test which is no longer appropriate.
tupek2 Oct 3, 2024
2126672
Apply style updates
Oct 3, 2024
06d2f03
Significantly simplify dual adjoint calculations.
tupek2 Oct 7, 2024
55c2b8a
Simplify reaction adjoint implementation even more.
tupek2 Oct 8, 2024
ab4ea41
Even more simplification of reaction adjoint loads, rename to setDual…
tupek2 Oct 8, 2024
d97afa9
Addition cleanup, reset sparse matrices to avoid highwater memory.
tupek2 Oct 8, 2024
a43765c
Test that reaction forces in qois are working when nonlinear and tran…
tupek2 Oct 8, 2024
396d791
Apply style updates
Oct 8, 2024
8efdcec
Small naming fix.
tupek2 Oct 8, 2024
6ed26e9
Fix docs.
tupek2 Oct 8, 2024
1a67141
Apply style updates
Oct 8, 2024
a8beeb6
Merge branch 'develop' into tupek/reaction_sensitivities
tupek2 Oct 9, 2024
9e6fc06
Merge branch 'develop' into tupek/reaction_sensitivities
tupek2 Oct 18, 2024
72c4868
Merge branch 'develop' into tupek/reaction_sensitivities
tupek2 Oct 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 27 additions & 55 deletions src/serac/physics/base_physics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_));
Expand Down Expand Up @@ -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
*
Expand Down Expand Up @@ -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<std::string, const serac::FiniteElementState&>)
{
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
Expand All @@ -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
Expand Down Expand Up @@ -550,7 +522,7 @@ class BasePhysics {
/**
* @brief List of adjoint finite element duals associated with this physics module
*/
std::vector<const serac::FiniteElementDual*> dual_adjoints_;
std::vector<const serac::FiniteElementState*> dual_adjoints_;

/// @brief The information needed for the physics parameters stored as Finite Element State fields
struct ParameterInfo {
Expand Down
139 changes: 75 additions & 64 deletions src/serac/physics/solid_mechanics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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<order, dim>{}, 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_,
Expand Down Expand Up @@ -216,7 +216,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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();
Expand Down Expand Up @@ -385,7 +385,8 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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;
Expand Down Expand Up @@ -753,15 +754,40 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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_;
}

/**
Expand Down Expand Up @@ -1134,8 +1160,10 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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_;
});
}
Expand Down Expand Up @@ -1201,7 +1229,9 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
std::unique_ptr<mfem::HypreParMatrix> m_mat(assemble(M));

// J = M + c0 * K
J_.reset();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this fixing a bug? Because the heat transfer module still does it the old way. If so, can you open an issue?

J_.reset(mfem::Add(1.0, *m_mat, c0_, *k_mat));
J_e_.reset();
J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_);

return *J_;
Expand Down Expand Up @@ -1316,9 +1346,11 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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()) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This check and the SLIC_ERROR_ROOT_IF above look like they're guarding the same problem. Can we combine these so that they don't get out of sync? Doesn't need to happen in this PR.

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");

Expand All @@ -1337,15 +1369,24 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
}
}

virtual void setDualAdjointBcs(std::unordered_map<std::string, const serac::FiniteElementState&> 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");
Expand All @@ -1361,11 +1402,20 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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<mfem::HypreParMatrix>(jacobian->Transpose());
J_.reset();
J_ = assemble(drdu);

auto J_T = std::unique_ptr<mfem::HypreParMatrix>(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);
Expand Down Expand Up @@ -1396,7 +1446,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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_;
Expand Down Expand Up @@ -1462,47 +1512,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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<mfem::HypreParMatrix> 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<DERIVATIVE>(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<DERIVATIVE>((*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<order, dim>;
Expand Down Expand Up @@ -1546,7 +1555,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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<ShapeAwareFunctional<shape_trial, test(trial, trial, parameter_space...)>> residual_;
Expand Down Expand Up @@ -1748,8 +1757,10 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, 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;

Expand Down
Loading