From eb65df3df76e901a7f1cd1b4f84486016be9e541 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Sat, 20 Jul 2024 06:59:17 -0700 Subject: [PATCH 1/9] Fix issue in warm start. For linear problems with body or traction loads, the warm start was solving the wrong linear system to start. Then we were solving another linear system again for no reason in the nonlinear solver after, greatly increasing costs in some cases. --- src/serac/physics/solid_mechanics.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 33ea4ab16..e55f33aa7 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1755,8 +1755,8 @@ class SolidMechanics, std::integer_se du_[j] -= displacement_(j); } - dr_ = 0.0; - mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, dr_); + dr_ = r; + dr_ *= -1.0; // Update the initial guess for changes in the parameters if this is not the first solve for (std::size_t parameter_index = 0; parameter_index < parameters_.size(); ++parameter_index) { @@ -1776,6 +1776,7 @@ class SolidMechanics, std::integer_se *parameters_[parameter_index].previous_state = *parameters_[parameter_index].state; } + mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, dr_); for (int i = 0; i < constrained_dofs.Size(); i++) { int j = constrained_dofs[i]; dr_[j] = du_[j]; @@ -1787,6 +1788,7 @@ class SolidMechanics, std::integer_se lin_solver.Mult(dr_, du_); displacement_ += du_; + } }; From 31b2b4a558065c75449fb4c91152313afa25696a Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Sat, 20 Jul 2024 15:53:59 -0700 Subject: [PATCH 2/9] style. --- src/serac/physics/solid_mechanics.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index e55f33aa7..faed9f458 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1788,7 +1788,6 @@ class SolidMechanics, std::integer_se lin_solver.Mult(dr_, du_); displacement_ += du_; - } }; From 3a856d4ee3d8821b60ca84a5ad5d34df798e0cee Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 1 Aug 2024 18:21:48 -0700 Subject: [PATCH 3/9] Simplify warm start, arguably make it more accurate. --- src/serac/physics/solid_mechanics.hpp | 57 ++++++++++++--------------- 1 file changed, 25 insertions(+), 32 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index faed9f458..8166deae0 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1641,19 +1641,6 @@ class SolidMechanics, std::integer_se displacement_, acceleration_, *parameters_[parameter_indices].state...); }...}; - /// @brief Array functions computing the derivative of the residual with respect to each given parameter evaluated at - /// the previous value of state - /// @note This is needed so the user can ask for a specific sensitivity at runtime as opposed to it being a - /// template parameter. - std::array< - std::function{}, 0.0, shape_displacement_, displacement_, acceleration_, - *parameters_[parameter_indices].previous_state...))(double)>, - sizeof...(parameter_indices)> - d_residual_d_previous_ = {[&](double t) { - return (*residual_)(DifferentiateWRT{}, t, shape_displacement_, - displacement_, acceleration_, *parameters_[parameter_indices].previous_state...); - }...}; - /// @brief Solve the Quasi-static Newton system virtual void quasiStaticSolve(double dt) { @@ -1734,12 +1721,36 @@ class SolidMechanics, std::integer_se /** * @brief Sets the Dirichlet BCs for the current time and computes an initial guess for parameters and displacement + * + * @note + * We want to solve + *\f$ + *r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) = 0 + *\f$ + *for $u_{n+1}$, given new values of parameters, essential b.c.s and time. The problem is that if we use $u_n$ as the initial guess for this new solve, most nonlinear solver algorithms will start off by linearizing at (or near) the initial guess. But, if the essential boundary conditions change by an amount on the order of the mesh size, then it's possible to invert elements and make that linearization point inadmissible (either because it converges slowly or that the inverted elements crash the program). + *So, we need a better initial guess. This "warm start" generates a guess by linear extrapolation from the previous known solution: + + *\f$ + *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx {r(u_n, p_n, U_n, t_n)} + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t + *\f$ + *If we assume that suddenly changing p and t will not lead to inverted elements, we can simplify the approximation to + *\f$ + *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dU} \Delta U + *\f$ + *Move all the known terms to the rhs and solve for $\Delta u$, + *\f$ + *\Delta u = - \bigg( \frac{dr_n}{du} \bigg)^{-1} \bigg( r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{dU} \Delta U \bigg) + *\f$ + *It is especially important to use the previously solved Jacobian in problems with material instabilities, as good nonlinear solvers will ensure positive definiteness at equilibrium. + *Once any parameter is changed, it is no longer certain to be positive definite, which will cause issues for many types linear solvers. */ void warmStartDisplacement(double dt) { // Update the linearized Jacobian matrix - auto [r, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].previous_state...); + auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, + *parameters_[parameter_indices].state...); J_ = assemble(drdu); J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); @@ -1758,24 +1769,6 @@ class SolidMechanics, std::integer_se dr_ = r; dr_ *= -1.0; - // Update the initial guess for changes in the parameters if this is not the first solve - for (std::size_t parameter_index = 0; parameter_index < parameters_.size(); ++parameter_index) { - // Compute the change in parameters parameter_diff = parameter_new - parameter_old - serac::FiniteElementState parameter_difference = *parameters_[parameter_index].state; - parameter_difference -= *parameters_[parameter_index].previous_state; - - // Compute a linearized estimate of the residual forces due to this change in parameter - auto drdparam = serac::get(d_residual_d_previous_[parameter_index](time_)); - auto residual_update = drdparam(parameter_difference); - - // Flip the sign to get the RHS of the Newton update system - // J^-1 du = - residual - dr_ -= residual_update; - - // Save the current parameter value for the next timestep - *parameters_[parameter_index].previous_state = *parameters_[parameter_index].state; - } - mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, dr_); for (int i = 0; i < constrained_dofs.Size(); i++) { int j = constrained_dofs[i]; From ef830f3f2eda28789bb1838fe9776b94de9d4995 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Fri, 2 Aug 2024 10:57:31 -0700 Subject: [PATCH 4/9] Get rid of some memory overhead \style. --- src/serac/physics/solid_mechanics.hpp | 14 ++++---------- src/serac/physics/solid_mechanics_contact.hpp | 1 - 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 8166deae0..b958f633b 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -259,7 +259,6 @@ class SolidMechanics, std::integer_se v_.SetSize(true_size); du_.SetSize(true_size); - dr_.SetSize(true_size); predicted_displacement_.SetSize(true_size); shape_displacement_ = 0.0; @@ -386,7 +385,6 @@ class SolidMechanics, std::integer_se u_ = 0.0; v_ = 0.0; du_ = 0.0; - dr_ = 0.0; predicted_displacement_ = 0.0; if (checkpoint_to_disk_) { @@ -1606,9 +1604,6 @@ class SolidMechanics, std::integer_se /// vector used to store the change in essential bcs between timesteps mfem::Vector du_; - /// vector used to store forces arising from du_ when applying time-dependent bcs - mfem::Vector dr_; - /// @brief used to communicate the ODE solver's predicted displacement to the residual operator mfem::Vector u_; @@ -1766,20 +1761,19 @@ class SolidMechanics, std::integer_se du_[j] -= displacement_(j); } - dr_ = r; - dr_ *= -1.0; + r *= -1.0; - mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, dr_); + mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, r); for (int i = 0; i < constrained_dofs.Size(); i++) { int j = constrained_dofs[i]; - dr_[j] = du_[j]; + r[j] = du_[j]; } auto& lin_solver = nonlin_solver_->linearSolver(); lin_solver.SetOperator(*J_); - lin_solver.Mult(dr_, du_); + lin_solver.Mult(r, du_); displacement_ += du_; } }; diff --git a/src/serac/physics/solid_mechanics_contact.hpp b/src/serac/physics/solid_mechanics_contact.hpp index 661ebed0e..17b80cf9b 100644 --- a/src/serac/physics/solid_mechanics_contact.hpp +++ b/src/serac/physics/solid_mechanics_contact.hpp @@ -253,7 +253,6 @@ class SolidMechanicsContact, using SolidMechanicsBase::d_residual_d_; using SolidMechanicsBase::DERIVATIVE; using SolidMechanicsBase::displacement_; - using SolidMechanicsBase::dr_; using SolidMechanicsBase::du_; using SolidMechanicsBase::J_; using SolidMechanicsBase::J_e_; From fef645590320017d26358f3136b18b223e550dac Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Fri, 2 Aug 2024 15:37:21 -0700 Subject: [PATCH 5/9] Fix style. --- src/serac/physics/solid_mechanics.hpp | 30 +++++++++++++++++---------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index b958f633b..d6c3bba6d 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1716,35 +1716,43 @@ class SolidMechanics, std::integer_se /** * @brief Sets the Dirichlet BCs for the current time and computes an initial guess for parameters and displacement - * + * * @note * We want to solve *\f$ *r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) = 0 *\f$ - *for $u_{n+1}$, given new values of parameters, essential b.c.s and time. The problem is that if we use $u_n$ as the initial guess for this new solve, most nonlinear solver algorithms will start off by linearizing at (or near) the initial guess. But, if the essential boundary conditions change by an amount on the order of the mesh size, then it's possible to invert elements and make that linearization point inadmissible (either because it converges slowly or that the inverted elements crash the program). - *So, we need a better initial guess. This "warm start" generates a guess by linear extrapolation from the previous known solution: + *for $u_{n+1}$, given new values of parameters, essential b.c.s and time. The problem is that if we use $u_n$ as the + initial guess for this new solve, most nonlinear solver algorithms will start off by linearizing at (or near) the + initial guess. But, if the essential boundary conditions change by an amount on the order of the mesh size, then it's + possible to invert elements and make that linearization point inadmissible (either because it converges slowly or + that the inverted elements crash the program). *So, we need a better initial guess. This "warm start" generates a + guess by linear extrapolation from the previous known solution: *\f$ - *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx {r(u_n, p_n, U_n, t_n)} + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t + *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx {r(u_n, p_n, U_n, t_n)} + \frac{dr_n}{du} \Delta u + + \frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t *\f$ *If we assume that suddenly changing p and t will not lead to inverted elements, we can simplify the approximation to *\f$ - *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dU} \Delta U + *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{du} \Delta u + + \frac{dr_n}{dU} \Delta U *\f$ *Move all the known terms to the rhs and solve for $\Delta u$, *\f$ - *\Delta u = - \bigg( \frac{dr_n}{du} \bigg)^{-1} \bigg( r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{dU} \Delta U \bigg) + *\Delta u = - \bigg( \frac{dr_n}{du} \bigg)^{-1} \bigg( r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{dU} \Delta U + \bigg) *\f$ - *It is especially important to use the previously solved Jacobian in problems with material instabilities, as good nonlinear solvers will ensure positive definiteness at equilibrium. - *Once any parameter is changed, it is no longer certain to be positive definite, which will cause issues for many types linear solvers. + *It is especially important to use the previously solved Jacobian in problems with material instabilities, as good + nonlinear solvers will ensure positive definiteness at equilibrium. *Once any parameter is changed, it is no longer + certain to be positive definite, which will cause issues for many types linear solvers. */ void warmStartDisplacement(double dt) { // Update the linearized Jacobian matrix auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].previous_state...); - auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, + auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, *parameters_[parameter_indices].state...); J_ = assemble(drdu); J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); @@ -1765,8 +1773,8 @@ class SolidMechanics, std::integer_se mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, r); for (int i = 0; i < constrained_dofs.Size(); i++) { - int j = constrained_dofs[i]; - r[j] = du_[j]; + int j = constrained_dofs[i]; + r[j] = du_[j]; } auto& lin_solver = nonlin_solver_->linearSolver(); From a37a2244268f994c411fe1a29b7b3dfdbcede223 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Fri, 2 Aug 2024 15:41:28 -0700 Subject: [PATCH 6/9] Fix latex. --- src/serac/physics/solid_mechanics.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index d6c3bba6d..592bd4e18 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1738,7 +1738,7 @@ class SolidMechanics, std::integer_se *0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dU} \Delta U *\f$ - *Move all the known terms to the rhs and solve for $\Delta u$, + *Move all the known terms to the rhs and solve for \f$\Delta u\f$, *\f$ *\Delta u = - \bigg( \frac{dr_n}{du} \bigg)^{-1} \bigg( r(u_n, p_{n+1}, U_n, t_{n+1}) + \frac{dr_n}{dU} \Delta U \bigg) From e9fe6f444621565cc402d135b47f108588ed3990 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 5 Aug 2024 12:34:52 -0700 Subject: [PATCH 7/9] Add an option to turn off warm start in solid mechanics. --- src/serac/physics/solid_mechanics.hpp | 56 +++++++++++++++------------ 1 file changed, 32 insertions(+), 24 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 592bd4e18..662e5f1aa 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -162,7 +162,7 @@ class SolidMechanics, std::integer_se SolidMechanics(std::unique_ptr solver, const serac::TimesteppingOptions timestepping_opts, const GeometricNonlinearities geom_nonlin, const std::string& physics_name, std::string mesh_tag, std::vector parameter_names = {}, int cycle = 0, double time = 0.0, - bool checkpoint_to_disk = false) + bool checkpoint_to_disk = false, bool use_warm_start = true) : BasePhysics(physics_name, mesh_tag, cycle, time, checkpoint_to_disk), displacement_( StateManager::newState(H1{}, detail::addPrefix(physics_name, "displacement"), mesh_tag_)), @@ -182,7 +182,8 @@ class SolidMechanics, std::integer_se ode2_(displacement_.space().TrueVSize(), {.time = ode_time_point_, .c0 = c0_, .c1 = c1_, .u = u_, .du_dt = v_, .d2u_dt2 = acceleration_}, *nonlin_solver_, bcs_), - geom_nonlin_(geom_nonlin) + geom_nonlin_(geom_nonlin), + use_warm_start_(use_warm_start) { SLIC_ERROR_ROOT_IF(mesh_.Dimension() != dim, axom::fmt::format("Compile time dimension, {0}, and runtime mesh dimension, {1}, mismatch", dim, @@ -1619,6 +1620,9 @@ class SolidMechanics, std::integer_se /// @brief A flag denoting whether to compute geometric nonlinearities in the residual GeometricNonlinearities geom_nonlin_; + /// @brief A flag denoting whether to compute the warm start for improved robustness + bool use_warm_start_; + /// @brief Coefficient containing the essential boundary values std::shared_ptr disp_bdr_coef_; @@ -1749,39 +1753,43 @@ class SolidMechanics, std::integer_se */ void warmStartDisplacement(double dt) { - // Update the linearized Jacobian matrix - auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].previous_state...); - auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, - *parameters_[parameter_indices].state...); - J_ = assemble(drdu); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); - du_ = 0.0; for (auto& bc : bcs_.essentials()) { // apply the future boundary conditions, but use the most recent Jacobians stiffness. bc.setDofs(du_, time_ + dt); } - auto& constrained_dofs = bcs_.allEssentialTrueDofs(); - for (int i = 0; i < constrained_dofs.Size(); i++) { - int j = constrained_dofs[i]; - du_[j] -= displacement_(j); - } + if (use_warm_start_) { + // Update the linearized Jacobian matrix + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, + *parameters_[parameter_indices].previous_state...); - r *= -1.0; + auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, + *parameters_[parameter_indices].state...); + J_ = assemble(drdu); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); - mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, r); - for (int i = 0; i < constrained_dofs.Size(); i++) { - int j = constrained_dofs[i]; - r[j] = du_[j]; - } + auto& constrained_dofs = bcs_.allEssentialTrueDofs(); + for (int i = 0; i < constrained_dofs.Size(); i++) { + int j = constrained_dofs[i]; + du_[j] -= displacement_(j); + } - auto& lin_solver = nonlin_solver_->linearSolver(); + r *= -1.0; - lin_solver.SetOperator(*J_); + mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du_, r); + for (int i = 0; i < constrained_dofs.Size(); i++) { + int j = constrained_dofs[i]; + r[j] = du_[j]; + } + + auto& lin_solver = nonlin_solver_->linearSolver(); + + lin_solver.SetOperator(*J_); + + lin_solver.Mult(r, du_); + } - lin_solver.Mult(r, du_); displacement_ += du_; } }; From 84626f1359f99b2a59c52bb521a9bc5354621b0a Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 5 Aug 2024 12:53:05 -0700 Subject: [PATCH 8/9] Try to reuse the jacobian when it exits. --- src/serac/physics/solid_mechanics.hpp | 29 +++++++++++++++++---------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 662e5f1aa..293b44a79 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -154,7 +154,8 @@ class SolidMechanics, std::integer_se * @param cycle The simulation cycle (i.e. timestep iteration) to intialize the physics module to * @param time The simulation time to initialize the physics module to * @param checkpoint_to_disk A flag to save the transient states on disk instead of memory for the transient adjoint - * solves + * @param use_warm_start A flag to turn on or off the displacement warm start predictor which helps robustness for + * large deformation problems * * @note On parallel file systems (e.g. lustre), significant slowdowns and occasional errors were observed when * writing and reading the needed trainsient states to disk for adjoint solves @@ -388,6 +389,9 @@ class SolidMechanics, std::integer_se du_ = 0.0; predicted_displacement_ = 0.0; + J_.reset(); + J_e_.reset(); + if (checkpoint_to_disk_) { outputStateToDisk(); } else { @@ -1759,20 +1763,23 @@ class SolidMechanics, std::integer_se bc.setDofs(du_, time_ + dt); } - if (use_warm_start_) { - // Update the linearized Jacobian matrix - auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].previous_state...); + auto& constrained_dofs = bcs_.allEssentialTrueDofs(); + for (int i = 0; i < constrained_dofs.Size(); i++) { + int j = constrained_dofs[i]; + du_[j] -= displacement_(j); + } + if (use_warm_start_ && is_quasistatic_) { + // Update the linearized Jacobian matrix auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, *parameters_[parameter_indices].state...); - J_ = assemble(drdu); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); - auto& constrained_dofs = bcs_.allEssentialTrueDofs(); - for (int i = 0; i < constrained_dofs.Size(); i++) { - int j = constrained_dofs[i]; - du_[j] -= displacement_(j); + // use the most recently evaluated Jacobian + if (!J_ || !J_e_) { + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, + *parameters_[parameter_indices].previous_state...); + J_ = assemble(drdu); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); } r *= -1.0; From b4eb80c9854240b9f03955ae9c63ff556bb59bda Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 5 Aug 2024 13:01:45 -0700 Subject: [PATCH 9/9] Always reassemble J for warm start for now. This can be optimized away later. --- src/serac/physics/solid_mechanics.hpp | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 293b44a79..5ac7c1663 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -210,9 +210,7 @@ class SolidMechanics, std::integer_se } adjoints_.push_back(&adjoint_displacement_); - duals_.push_back(&reactions_); - dual_adjoints_.push_back(&reactions_adjoint_load_); // Create a pack of the primal field and parameter finite element spaces @@ -259,7 +257,6 @@ class SolidMechanics, std::integer_se u_.SetSize(true_size); v_.SetSize(true_size); - du_.SetSize(true_size); predicted_displacement_.SetSize(true_size); @@ -389,9 +386,6 @@ class SolidMechanics, std::integer_se du_ = 0.0; predicted_displacement_ = 0.0; - J_.reset(); - J_e_.reset(); - if (checkpoint_to_disk_) { outputStateToDisk(); } else { @@ -1775,12 +1769,10 @@ class SolidMechanics, std::integer_se *parameters_[parameter_indices].state...); // use the most recently evaluated Jacobian - if (!J_ || !J_e_) { - auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].previous_state...); - J_ = assemble(drdu); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); - } + auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, + *parameters_[parameter_indices].previous_state...); + J_ = assemble(drdu); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); r *= -1.0;