diff --git a/ContinuationSolvers b/ContinuationSolvers index b7c26c473f..c59883d156 160000 --- a/ContinuationSolvers +++ b/ContinuationSolvers @@ -1 +1 @@ -Subproject commit b7c26c473ffb1aaaefc82832b6e1f625b239992c +Subproject commit c59883d15679cd764fe88865c06c41469854fc87 diff --git a/examples/contact/homotopy/two_blocks.cpp b/examples/contact/homotopy/two_blocks.cpp index a37df57c21..570456436b 100644 --- a/examples/contact/homotopy/two_blocks.cpp +++ b/examples/contact/homotopy/two_blocks.cpp @@ -73,10 +73,11 @@ class TiedContactProblem : public EqualityConstrainedHomotopyProblem { std::shared_ptr constraints, mfem::Array fixed_tdof_list, mfem::Array disp_tdof_list, mfem::Vector uDC); mfem::Vector residual(const mfem::Vector& u, bool fresh_evaluation) const; - mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool fresh_evaluation); + mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool fresh_evaluation, bool fresh_derivative); mfem::Vector constraint(const mfem::Vector& u, bool fresh_evaluation) const; - mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool fresh_evaluation); - mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool fresh_evaluation) const; + mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool fresh_evaluation, bool fresh_derivative); + mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool fresh_evaluation, + bool fresh_derivative) const; void fullDisplacement(const mfem::Vector& x, mfem::Vector& u); virtual ~TiedContactProblem(); }; @@ -290,7 +291,8 @@ int main(int argc, char* argv[]) SLIC_WARNING_ROOT_IF(!converged, "Homotopy solver did not converge"); // visualize - auto writer = createParaviewOutput(mesh->mfemParMesh(), smith::getConstFieldPointers(states), "contact"); + auto writer = + createParaviewOutput(mesh->mfemParMesh(), smith::getConstFieldPointers(states), "two_block_tiedcontact_plot"); if (visualize) { mfem::Vector u(states[FIELD::DISP].space().GetTrueVSize()); u = problem.GetDisplacement(X0); @@ -349,6 +351,7 @@ TiedContactProblem::TiedContactProblem(std::vector mfem::Vector TiedContactProblem::residual(const mfem::Vector& u, bool /*fresh_evaluation*/) const { + // TODO: cache this evaluation residual_states_[FIELD::DISP]->Set(1.0, u); auto res = weak_form_->residual(time_info_, shape_disp_.get(), smith::getConstFieldPointers(residual_states_)); return res; @@ -356,18 +359,23 @@ mfem::Vector TiedContactProblem::residual(const mfem::Vector& template mfem::HypreParMatrix* TiedContactProblem::residualJacobian(const mfem::Vector& u, - bool /*fresh_evaluation*/) + bool fresh_evaluation, + bool fresh_derivative) { - residual_states_[FIELD::DISP]->Set(1.0, u); - drdu_ = weak_form_->jacobian(time_info_, shape_disp_.get(), smith::getConstFieldPointers(residual_states_), - jacobian_weights_); + if (fresh_evaluation || fresh_derivative) { + residual_states_[FIELD::DISP]->Set(1.0, u); + drdu_ = weak_form_->jacobian(time_info_, shape_disp_.get(), smith::getConstFieldPointers(residual_states_), + jacobian_weights_); + } else { + SLIC_ERROR_ROOT_IF(!drdu_, "must call evaluation method prior to using cached data"); + } return drdu_.get(); } template -mfem::Vector TiedContactProblem::constraint(const mfem::Vector& u, bool /*fresh_evaluation*/) const +mfem::Vector TiedContactProblem::constraint(const mfem::Vector& u, bool fresh_evaluation) const { - bool fresh_evaluation = true; + // No caching necessary as that is already done in ContactConstraint::evaluate contact_states_[smith::ContactFields::DISP]->Set(1.0, u); auto gap = constraints_->evaluate(time_info_.time(), time_info_.dt(), smith::getConstFieldPointers(contact_states_), fresh_evaluation); @@ -376,24 +384,29 @@ mfem::Vector TiedContactProblem::constraint(const mfem::Vecto template mfem::HypreParMatrix* TiedContactProblem::constraintJacobian(const mfem::Vector& u, - bool /*fresh_evaluation*/) + bool fresh_evaluation, + bool fresh_derivative) { - bool fresh_evaluation = true; - contact_states_[smith::ContactFields::DISP]->Set(1.0, u); - dcdu_ = constraints_->jacobian(time_info_.time(), time_info_.dt(), smith::getConstFieldPointers(contact_states_), - smith::ContactFields::DISP, fresh_evaluation); + if (fresh_evaluation || fresh_derivative) { + contact_states_[smith::ContactFields::DISP]->Set(1.0, u); + dcdu_ = constraints_->jacobian(time_info_.time(), time_info_.dt(), smith::getConstFieldPointers(contact_states_), + smith::ContactFields::DISP, fresh_evaluation, fresh_derivative); + } else { + SLIC_ERROR_ROOT_IF(!dcdu_, "must call evaluation method prior to using cached data"); + } return dcdu_.get(); } template mfem::Vector TiedContactProblem::constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, - bool /*fresh_evaluation*/) const + bool fresh_evaluation, + bool fresh_derivative) const { - bool fresh_evaluation = true; + // internal caching done in ContactConstraint::residual_contribution contact_states_[smith::ContactFields::DISP]->Set(1.0, u); - auto res_contribution = constraints_->residual_contribution(time_info_.time(), time_info_.dt(), - smith::getConstFieldPointers(contact_states_), l, - smith::ContactFields::DISP, fresh_evaluation); + auto res_contribution = constraints_->residual_contribution( + time_info_.time(), time_info_.dt(), smith::getConstFieldPointers(contact_states_), l, smith::ContactFields::DISP, + fresh_evaluation, fresh_derivative); return res_contribution; } diff --git a/examples/inertia_relief/inertia_relief_example.cpp b/examples/inertia_relief/inertia_relief_example.cpp index 5aaba18fe3..9744e991b4 100644 --- a/examples/inertia_relief/inertia_relief_example.cpp +++ b/examples/inertia_relief/inertia_relief_example.cpp @@ -117,10 +117,11 @@ class InertialReliefProblem : public EqualityConstrainedHomotopyProblem { std::shared_ptr mesh, std::shared_ptr weak_form, std::vector> constraints); mfem::Vector residual(const mfem::Vector& u, bool fresh_evaluation) const; - mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool fresh_evaluation) const; + mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool fresh_evaluation, + bool fresh_derivative) const; mfem::Vector constraint(const mfem::Vector& u, bool fresh_evaluation) const; - mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool fresh_evaluation); - mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool fresh_evaluation); + mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool fresh_evaluation, bool fresh_derivative); + mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool fresh_evaluation, bool fresh_derivative); virtual ~InertialReliefProblem(); }; @@ -362,9 +363,9 @@ mfem::Vector InertialReliefProblem::residual(const mfem::Vector& u, bool fresh_e // constraint Jacobian transpose vector product mfem::Vector InertialReliefProblem::constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, - bool fresh_evaluation) const + bool fresh_evaluation, bool fresh_derivative) const { - if (fresh_evaluation) { + if (fresh_evaluation || fresh_derivative) { int dim_constraints = GetMultiplierDim(); int dim_displacement = GetDisplacementDim(); obj_states_[FIELD::DISP]->Set(1.0, u); @@ -389,9 +390,10 @@ mfem::Vector InertialReliefProblem::constraintJacobianTvp(const mfem::Vector& u, } // Jacobian of the residual -mfem::HypreParMatrix* InertialReliefProblem::residualJacobian(const mfem::Vector& u, bool fresh_evaluation) +mfem::HypreParMatrix* InertialReliefProblem::residualJacobian(const mfem::Vector& u, bool fresh_evaluation, + bool fresh_derivative) { - if (fresh_evaluation) { + if (fresh_evaluation || fresh_derivative) { obj_states_[FIELD::DISP]->Set(1.0, u); drdu_.reset(); drdu_ = weak_form_->jacobian(time_info_, shape_disp_.get(), getConstFieldPointers(all_states_), jacobian_weights_); @@ -425,11 +427,12 @@ mfem::Vector InertialReliefProblem::constraint(const mfem::Vector& u, bool fresh } // Jacobian of the constraint -mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vector& u, bool fresh_evaluation) +mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vector& u, bool fresh_evaluation, + bool fresh_derivative) { - int dim_constraints = GetMultiplierDim(); - int glbdim_displacement = GetGlobalDisplacementDim(); - if (fresh_evaluation) { + if (fresh_evaluation || fresh_derivative) { + int dim_constraints = GetMultiplierDim(); + int glbdim_displacement = GetGlobalDisplacementDim(); obj_states_[FIELD::DISP]->Set(1.0, u); // dense rows int nentries = glbdim_displacement; diff --git a/src/smith/physics/constraint.hpp b/src/smith/physics/constraint.hpp index 6d58ac8790..20335cd471 100644 --- a/src/smith/physics/constraint.hpp +++ b/src/smith/physics/constraint.hpp @@ -54,11 +54,14 @@ class Constraint { * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::unique_ptr */ virtual std::unique_ptr jacobian(double time, double dt, const std::vector& fields, int direction, - bool fresh_evaluation = true) const = 0; + bool fresh_evaluation = true, + bool fresh_derivative = true) const = 0; /** @brief Virtual interface for computing constraint Jacobian_tilde from a vector of smith::FiniteElementState* * Jacobian_tilde is an optional approximation of the true Jacobian @@ -68,13 +71,16 @@ class Constraint { * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::unique_ptr */ virtual std::unique_ptr jacobian_tilde(double time, double dt, const std::vector& fields, int direction, - bool fresh_evaluation = true) const + bool fresh_evaluation = true, + bool fresh_derivative = true) const { - return jacobian(time, dt, fields, direction, fresh_evaluation); + return jacobian(time, dt, fields, direction, fresh_evaluation, fresh_derivative); }; /** @brief Virtual interface for computing residual contribution Jacobian_tilde^(Transpose) * (Lagrange multiplier) @@ -86,13 +92,16 @@ class Constraint { * @param multipliers mfem::Vector of Lagrange multipliers * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::Vector */ virtual mfem::Vector residual_contribution(double time, double dt, const std::vector& fields, const mfem::Vector& multipliers, int direction, - bool fresh_evaluation = true) const + bool fresh_evaluation = true, bool fresh_derivative = true) const { - std::unique_ptr jac = jacobian_tilde(time, dt, fields, direction, fresh_evaluation); + std::unique_ptr jac = + jacobian_tilde(time, dt, fields, direction, fresh_evaluation, fresh_derivative); mfem::Vector y(jac->Width()); y = 0.0; SLIC_ERROR_ROOT_IF(jac->Height() != multipliers.Size(), "Incompatible matrix and vector sizes."); @@ -109,12 +118,15 @@ class Constraint { * @param multipliers mfem::Vector of Lagrange multipliers * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::unique_ptr */ virtual std::unique_ptr residual_contribution_jacobian( [[maybe_unused]] double time, [[maybe_unused]] double dt, [[maybe_unused]] const std::vector& fields, [[maybe_unused]] const mfem::Vector& multipliers, - [[maybe_unused]] int direction, [[maybe_unused]] bool fresh_evaluation = true) const + [[maybe_unused]] int direction, [[maybe_unused]] bool fresh_evaluation = true, + [[maybe_unused]] bool fresh_derivative = true) const { SLIC_ERROR_ROOT(axom::fmt::format("Base class must override residual_contribution_jacobian before usage")); std::unique_ptr res_contr_jacobian = nullptr; diff --git a/src/smith/physics/contact_constraint.hpp b/src/smith/physics/contact_constraint.hpp index 4045aa0079..6d0b2439c7 100644 --- a/src/smith/physics/contact_constraint.hpp +++ b/src/smith/physics/contact_constraint.hpp @@ -44,41 +44,23 @@ enum ContactFields /** @brief Interface for extracting the iblock, jblock block from a std::unique_ptr * said block is returned as a std::unique_ptr if possible - * All other blocks are deleted * * @param block_operator the block operator * @param iblock row block index * @param jblock column block index * @return std::unique_ptr The requested block of block_operator */ -static std::unique_ptr safelyObtainBlock(mfem::BlockOperator* block_operator, int iblock, - int jblock, bool own_blocks = false) +static std::unique_ptr obtainBlock(mfem::BlockOperator* block_operator, int iblock, int jblock) { SLIC_ERROR_IF(iblock < 0 || jblock < 0, "block indicies must be non-negative"); SLIC_ERROR_IF(iblock > block_operator->NumRowBlocks() || jblock > block_operator->NumColBlocks(), "one or more block indicies are too large and the requested block does not exist"); SLIC_ERROR_IF(block_operator->IsZeroBlock(iblock, jblock), "attempting to extract a null block"); - block_operator->owns_blocks = false; - for (int i = 0; i < block_operator->NumRowBlocks(); i++) { - for (int j = 0; j < block_operator->NumColBlocks(); j++) { - if (i == iblock && j == jblock) { - continue; - } - if (!block_operator->IsZeroBlock(i, j) && !own_blocks) { - delete &block_operator->GetBlock(i, j); - } - } - } auto Ablk = dynamic_cast(&block_operator->GetBlock(iblock, jblock)); SLIC_ERROR_IF(!Ablk, "failed cast block to mfem::HypreParMatrix"); - if (own_blocks) { - // deep copy --> unique_ptr - auto Ablk_unique = std::make_unique(*Ablk); - return Ablk_unique; - } else { - std::unique_ptr Ablk_unique(Ablk); - return Ablk_unique; - } + // deep copy --> unique_ptr + auto Ablk_unique = std::make_unique(*Ablk); + return Ablk_unique; }; class FiniteElementState; @@ -111,6 +93,7 @@ class ContactConstraint : public Constraint { contact_opts_.enforcement = ContactEnforcement::LagrangeMultiplier; contact_.addContactInteraction(interaction_id, bdry_attr_surf1, bdry_attr_surf2, contact_opts_); interaction_id_ = interaction_id; + tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP); } /// @brief destructor @@ -128,20 +111,18 @@ class ContactConstraint : public Constraint { mfem::Vector evaluate(double time, double dt, const std::vector& fields, bool fresh_evaluation = true) const override { - contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); - tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP); - if (fresh_evaluation) { + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); // note: Tribol does not use cycle. int cycle = 0; contact_.update(cycle, time, dt); - pressures_set_ = false; + auto gaps_hpv = contact_.mergedGaps(false); + // Note: this copy is needed to prevent the HypreParVector pointer from going out of scope. see + // https://github.com/mfem/mfem/issues/5029 + cached_gap_.SetSize(gaps_hpv.Size()); + cached_gap_.Set(1.0, gaps_hpv); } - auto gaps_hpv = contact_.mergedGaps(false); - // Note: this copy is needed to prevent the HypreParVector pointer from going out of scope. see - // https://github.com/mfem/mfem/issues/5029 - mfem::Vector gaps = gaps_hpv; - return gaps; + return cached_gap_; }; /** @brief Interface for computing contact gap constraint Jacobian from a vector of smith::FiniteElementState* @@ -151,24 +132,27 @@ class ContactConstraint : public Constraint { * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::unique_ptr The true Jacobian */ std::unique_ptr jacobian(double time, double dt, const std::vector& fields, - int direction, bool fresh_evaluation = true) const override + int direction, bool fresh_evaluation = true, + [[maybe_unused]] bool fresh_derivative = true) const override { SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); - contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); - tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); - + /* If this is not a new evaluation point e.g., + if the constraint was evaluated at this point + then the derivatives have already been computed + and we do not need to recompute them.*/ if (fresh_evaluation) { + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); int cycle = 0; contact_.update(cycle, time, dt); J_contact_ = contact_.mergedJacobian(); - pressures_set_ = false; } - int iblock = 1; - int jblock = 0; - auto dgdu = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + // obtain (1, 0) block entry from the 2 x 2 block contact linear system + auto dgdu = obtainBlock(J_contact_.get(), 1, 0); return dgdu; }; @@ -181,36 +165,33 @@ class ContactConstraint : public Constraint { * @param multipliers mfem::Vector of Lagrange multipliers * @param direction index for which field to take the gradient with respect to * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation + * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached + * evaluation * @return std::Vector */ mfem::Vector residual_contribution(double time, double dt, const std::vector& fields, - const mfem::Vector& multipliers, int direction, - bool fresh_evaluation = true) const override + const mfem::Vector& multipliers, int direction, bool fresh_evaluation = true, + bool fresh_derivative = true) const override { SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); - contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); - tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP); - int cycle = 0; if (fresh_evaluation) { - // we need to call update first to update gaps + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); + // first update gaps for (auto& interaction : contact_.getContactInteractions()) { interaction.evalJacobian(false); } contact_.update(cycle, time, dt); - pressures_set_ = false; } - if (!pressures_set_) { - // with updated gaps, we can update pressure for contact interactions with penalty enforcement + if (fresh_evaluation || fresh_derivative) { + // with updated gaps, then update pressure for contact interactions with penalty enforcement contact_.setPressures(multipliers); // call update again with the right pressures for (auto& interaction : contact_.getContactInteractions()) { interaction.evalJacobian(true); } contact_.update(cycle, time, dt); - pressures_set_ = true; } - return contact_.forces(); }; @@ -228,22 +209,21 @@ class ContactConstraint : public Constraint { std::unique_ptr residual_contribution_jacobian(double time, double dt, const std::vector& fields, const mfem::Vector& multipliers, int direction, - bool fresh_evaluation = true) const override + bool fresh_evaluation = true, + bool fresh_derivative = true) const override { SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); - contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); - tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); int cycle = 0; if (fresh_evaluation) { - // we need to call update first to update gaps + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); + // first update gaps for (auto& interaction : contact_.getContactInteractions()) { interaction.evalJacobian(false); } contact_.update(cycle, time, dt); - pressures_set_ = false; } - if (!pressures_set_) { + if (fresh_evaluation || fresh_derivative) { // with updated gaps, we can update pressure for contact interactions with penalty enforcement contact_.setPressures(multipliers); // call update again with the right pressures @@ -252,11 +232,9 @@ class ContactConstraint : public Constraint { } contact_.update(cycle, time, dt); J_contact_ = contact_.mergedJacobian(); - pressures_set_ = true; } - int iblock = 0; - int jblock = 0; - auto Hessian = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + // obtain (0, 0) block entry from the 2 x 2 block contact linear system + auto Hessian = obtainBlock(J_contact_.get(), 0, 0); return Hessian; }; @@ -270,22 +248,19 @@ class ContactConstraint : public Constraint { * @return std::unique_ptr */ std::unique_ptr jacobian_tilde(double time, double dt, const std::vector& fields, - int direction, bool fresh_evaluation = true) const override + int direction, bool fresh_evaluation = true, + [[maybe_unused]] bool fresh_derivative = true) const override { SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); - contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); - tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); int cycle = 0; if (fresh_evaluation) { + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); contact_.update(cycle, time, dt); - J_contact_.reset(); J_contact_ = contact_.mergedJacobian(); - pressures_set_ = false; } - int iblock = 0; - int jblock = 1; - auto dgduT = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + // obtain (0, 1) block entry from the 2 x 2 block contact linear system + auto dgduT = obtainBlock(J_contact_.get(), 0, 1); std::unique_ptr dgdu(dgduT->Transpose()); return dgdu; }; @@ -314,14 +289,9 @@ class ContactConstraint : public Constraint { mutable std::unique_ptr J_contact_; /** - * @brief own_blocks_ temporary boolean - */ - const bool own_blocks_ = true; - - /** - * @brief pressures_set_ are the Lagrange multipliers set + * @brief cached_gap_ gap function cached in order to not redo any unnecessary computations */ - mutable bool pressures_set_ = false; + mutable mfem::Vector cached_gap_; }; } // namespace smith