Put in a recursive strategy to slowly ramp up new boundary conditions…

src/serac/numerics/equation_solver.cpp

@@ -19,6 +19,20 @@
 
 namespace serac {
 
+/// @brief A simple exception class for recording Newton and TrustRegion failures
+class SolverException : public std::exception {
+ public:
+  /// constructor
+  SolverException(const std::string& message) : msg(message) {}
+
+  /// the message returned by the exception
+  const char* what() const noexcept override { return msg.c_str(); }
+
+ private:
+  /// message string
+  std::string msg;
+};
+
 /// Newton solver with a 2-way line-search.  Reverts to regular Newton if max_line_search_iterations is set to 0.
 class NewtonSolver : public mfem::NewtonSolver {
  protected:
@@ -108,8 +122,11 @@ class NewtonSolver : public mfem::NewtonSolver {
       }
 
       if (norm != norm) {
-        mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
-        mfem::out << "Newton: No convergence!\n";
+        if (print_options.first_and_last) {
+          mfem::out << "Initial residual for Newton iteration is undefined/nan." << std::endl;
+          mfem::out << "Newton: No convergence!\n";
+        }
+        throw SolverException("Initial residual for Newton iteration is undefined/nan.");
         return;
       }
 
@@ -199,6 +216,9 @@ class NewtonSolver : public mfem::NewtonSolver {
     if (!converged && (print_options.summary || print_options.warnings)) {
       mfem::out << "Newton: No convergence!\n";
     }
+    if (!converged) {
+      throw SolverException("Newton algorithm failed to converge.");
+    }
   }
 };
 
@@ -665,14 +685,17 @@ class TrustRegion : public mfem::NewtonSolver {
           mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm);
           mfem::out << ", x_incr = " << std::setw(13) << trResults.d.Norml2();
         } else {
-          mfem::out << ", norm goal = " << std::setw(13) << norm_goal << "\n";
+          mfem::out << ", norm goal = " << std::setw(13) << norm_goal;
         }
-        mfem::out << '\n';
+        mfem::out << "\n";
       }
 
       if (norm != norm) {
-        mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
-        mfem::out << "Newton: No convergence!\n";
+        if (print_options.first_and_last) {
+          mfem::out << "Initial residual for trust-region iteration is undefined/nan." << std::endl;
+          mfem::out << "Trust-Region: No convergence!\n";
+        }
+        throw SolverException("Initial residual for trust-region iteration is undefined/nan.");
         return;
       }
 
@@ -791,8 +814,8 @@ class TrustRegion : public mfem::NewtonSolver {
           happyAboutTrSize = true;
           if (print_options.iterations) {
             printTrustRegionInfo(realObjective, modelObjective, trResults.cg_iterations_count, tr_size, true);
-            trResults.cg_iterations_count =
-                0;  // zero this output so it doesn't look like the linesearch is doing cg iterations
+            // zero this output so it doesn't look like the linesearch is doing cg iterations
+            trResults.cg_iterations_count = 0;
           }
           break;
         }
@@ -863,6 +886,10 @@ class TrustRegion : public mfem::NewtonSolver {
       mfem::out << "num subspace solves = " << num_subspace_solves << "\n";
       mfem::out << "num jacobian_assembles = " << num_jacobian_assembles << "\n";
     }
+
+    if (!converged) {
+      throw SolverException("Trust-region algorithm failed to converge.");
+    }
   }
 };
 

src/serac/physics/solid_mechanics.hpp

@@ -131,17 +131,20 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    * @param checkpoint_to_disk Flag to save the transient states on disk instead of memory for transient adjoint solver
    * @param use_warm_start Flag to turn on or off the displacement warm start predictor which helps robustness for
    * large deformation problems
+   * @param max_timestep_cutbacks The maximum number of times the supplied timestep can be cutback before the solver
+   * gives up
    *
    * @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
    */
   SolidMechanics(const NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions lin_opts,
                  const serac::TimesteppingOptions timestepping_opts, const std::string& physics_name,
                  std::string mesh_tag, std::vector<std::string> parameter_names = {}, int cycle = 0, double time = 0.0,
-                 bool checkpoint_to_disk = false, bool use_warm_start = true)
+                 bool checkpoint_to_disk = false, bool use_warm_start = true, int max_timestep_cutbacks = 5)
       : SolidMechanics(
             std::make_unique<EquationSolver>(nonlinear_opts, lin_opts, StateManager::mesh(mesh_tag).GetComm()),
-            timestepping_opts, physics_name, mesh_tag, parameter_names, cycle, time, checkpoint_to_disk, use_warm_start)
+            timestepping_opts, physics_name, mesh_tag, parameter_names, cycle, time, checkpoint_to_disk, use_warm_start,
+            max_timestep_cutbacks)
   {
   }
 
@@ -158,13 +161,16 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    * @param checkpoint_to_disk Flag to save the transient states on disk instead of memory for 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
+   * @param max_timestep_cutbacks The maximum number of times the supplied timestep can be cutback before the solver
+   * gives up
    *
    * @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
    */
   SolidMechanics(std::unique_ptr<serac::EquationSolver> solver, const serac::TimesteppingOptions timestepping_opts,
                  const std::string& physics_name, std::string mesh_tag, std::vector<std::string> parameter_names = {},
-                 int cycle = 0, double time = 0.0, bool checkpoint_to_disk = false, bool use_warm_start = true)
+                 int cycle = 0, double time = 0.0, bool checkpoint_to_disk = false, bool use_warm_start = true,
+                 int max_timestep_cutbacks = 5)
       : BasePhysics(physics_name, mesh_tag, cycle, time, checkpoint_to_disk),
         displacement_(
             StateManager::newState(H1<order, dim>{}, detail::addPrefix(physics_name, "displacement"), mesh_tag_)),
@@ -184,7 +190,8 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
         ode2_(displacement_.space().TrueVSize(),
               {.time = time_, .c0 = c0_, .c1 = c1_, .u = u_, .du_dt = v_, .d2u_dt2 = acceleration_}, *nonlin_solver_,
               bcs_),
-        use_warm_start_(use_warm_start)
+        use_warm_start_(use_warm_start),
+        max_timestep_cutbacks_(max_timestep_cutbacks)
   {
     SERAC_MARK_FUNCTION;
     SLIC_ERROR_ROOT_IF(mesh_.Dimension() != dim,
@@ -1176,7 +1183,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
     }
 
     if (is_quasistatic_) {
-      quasiStaticSolve(dt);
+      quasiStaticSolve(dt, 1.0, 0);
     } else {
       // The current ode interface tracks 2 times, one internally which we have a handle to via time_,
       // and one here via the step interface.
@@ -1500,6 +1507,9 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   /// @brief A flag denoting whether to compute the warm start for improved robustness
   bool use_warm_start_;
 
+  /// @brief The maximum number of cut-back in the supplied time-step increment before the solver will give up
+  int max_timestep_cutbacks_;
+
   /// @brief Coefficient containing the essential boundary values
   std::shared_ptr<mfem::VectorCoefficient> disp_bdr_coef_;
 
@@ -1518,16 +1528,46 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       }...};
 
   /// @brief Solve the Quasi-static Newton system
-  virtual void quasiStaticSolve(double dt)
+  virtual void quasiStaticSolve(double dt, double step_fraction_of_dt_to_try, int level)
   {
-    // warm start must be called prior to the time update so that the previous Jacobians can be used consistently
-    // throughout.
-    warmStartDisplacement(dt);
-    time_ += dt;
-
-    // this method is essentially equivalent to the 1-liner
-    // u += dot(inv(J), dot(J_elim[:, dofs], (U(t + dt) - u)[dofs]));
-    nonlin_solver_->solve(displacement_);
+    if (level > max_timestep_cutbacks_) {
+      if (mpi_rank_ == 0)
+        std::cout << "Too many boundary condition cutbacks, accepting solution even though there may be issues. Try "
+                     "increasing the number of load steps "
+                  << std::endl;
+      return;
+    }
+
+    static constexpr bool solver_info_print = false;
+
+    bool solver_success = true;
+    try {
+      // warm start must be called prior to the time update so that the previous Jacobians can be used consistently
+      // throughout.
+      if (level == 0) time_ += dt;
+      warmStartDisplacement(dt, step_fraction_of_dt_to_try);
+      nonlin_solver_->solve(displacement_);
+    } catch (const std::exception& e) {
+      if (mpi_rank_ == 0 && solver_info_print) {
+        mfem::out << "Caught nonlinear solver exception: " << e.what() << std::endl;
+      }
+      displacement_ -= du_;
+      solver_success = false;
+      quasiStaticSolve(dt, 0.5 * step_fraction_of_dt_to_try, level + 1);
+      quasiStaticSolve(dt, 1.0, level + 1);
+    }
+
+    if (solver_success) {
+      if (step_fraction_of_dt_to_try == 1.0) {
+        if (mpi_rank_ == 0 && solver_info_print) {
+          mfem::out << "SolidMechanics solve succeeded for time " << time_ << " dt = " << dt << std::endl;
+        }
+      } else {
+        if (mpi_rank_ == 0 && solver_info_print) {
+          mfem::out << "SolidMechanics substep solve succeeded for time " << time_ << " dt = " << dt << std::endl;
+        }
+      }
+    }
   }
 
   /**
@@ -1625,14 +1665,14 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
    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)
+  void warmStartDisplacement(double dt, double displacement_scale_factor)
   {
     SERAC_MARK_FUNCTION;
 
     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);
+      bc.setDofs(du_, time_);
     }
 
     auto& constrained_dofs = bcs_.allEssentialTrueDofs();
@@ -1641,13 +1681,17 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       du_[j] -= displacement_(j);
     }
 
+    if (displacement_scale_factor != 1.0) {
+      du_ *= displacement_scale_factor;
+    }
+
     if (use_warm_start_) {
       // Update the linearized Jacobian matrix
-      auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_,
+      auto r = (*residual_)(time_, shape_displacement_, displacement_, acceleration_,
                             *parameters_[parameter_indices].state...);
 
       // use the most recently evaluated Jacobian
-      auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
+      auto [_, drdu] = (*residual_)(time_ - dt, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
                                     *parameters_[parameter_indices].previous_state...);
       J_.reset();
       J_ = assemble(drdu);
@@ -1665,11 +1709,23 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       auto& lin_solver = nonlin_solver_->linearSolver();
 
       lin_solver.SetOperator(*J_);
-
       lin_solver.Mult(r, du_);
     }
 
     displacement_ += du_;
+
+    // due to solver inaccuracies, the end of step boundary conditions may not be exact at this point
+    if (displacement_scale_factor == 1.0) {
+      du_ = 0.0;
+      for (auto& bc : bcs_.essentials()) {
+        bc.setDofs(du_, time_);
+      }
+      for (int i = 0; i < constrained_dofs.Size(); i++) {
+        int j = constrained_dofs[i];
+        du_[j] -= displacement_(j);
+      }
+      displacement_ += du_;
+    }
   }
 };
 

src/serac/physics/solid_mechanics_contact.hpp

@@ -221,7 +221,8 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
 
  protected:
   /// @brief Solve the Quasi-static Newton system
-  void quasiStaticSolve(double dt) override
+  void quasiStaticSolve(double dt, [[maybe_unused]] double step_fraction_of_dt_to_try,
+                        [[maybe_unused]] int level) override
   {
     // warm start must be called prior to the time update so that the previous Jacobians can be used consistently
     // throughout. warm start for contact needs to include the previous stiffness terms associated with contact

src/serac/physics/tests/contact_patch_tied.cpp

@@ -70,9 +70,9 @@ TEST_P(ContactPatchTied, patch)
 
   NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton,
                                            .relative_tol = 1.0e-10,
-                                           .absolute_tol = 1.0e-10,
-                                           .max_iterations = 20,
-                                           .print_level = 1};
+                                           .absolute_tol = 5.0e-10,
+                                           .max_iterations = 25,
+                                           .print_level = 0};
 
   ContactOptions contact_options{.method = ContactMethod::SingleMortar,
                                  .enforcement = GetParam().first,

src/serac/physics/tests/lce_Brighenti_tensile.cpp

@@ -80,7 +80,7 @@ TEST(LiquidCrystalElastomer, Brighenti)
 
   // Construct a solid mechanics solver
   LinearSolverOptions linear_options = {
-      .linear_solver = LinearSolver::GMRES,
+      .linear_solver = LinearSolver::CG,
       .preconditioner = Preconditioner::HypreAMG,
       .relative_tol = 1.0e-6,
       .absolute_tol = 1.0e-14,

src/serac/physics/tests/solid_statics_patch.cpp

@@ -307,22 +307,21 @@ double pressure_error()
   Domain material_block = EntireDomain(pmesh);
   Domain boundary = EntireBoundary(pmesh);
 
-// Construct a solid mechanics solver
-#ifdef SERAC_USE_SUNDIALS
-  serac::NonlinearSolverOptions nonlin_solver_options{.nonlin_solver = NonlinearSolver::KINBacktrackingLineSearch,
-                                                      .relative_tol = 0.0,
-                                                      .absolute_tol = 1.0e-14,
-                                                      .max_iterations = 30};
-#else
+  // Construct a solid mechanics solver
   serac::NonlinearSolverOptions nonlin_solver_options{
       .nonlin_solver = NonlinearSolver::Newton, .relative_tol = 0.0, .absolute_tol = 1.0e-14, .max_iterations = 30};
-#endif
 
-  auto equation_solver = std::make_unique<EquationSolver>(
-      nonlin_solver_options, serac::solid_mechanics::default_linear_options, pmesh.GetComm());
+  LinearSolverOptions cg_solver_options{.linear_solver = LinearSolver::CG,
+                                        .preconditioner = Preconditioner::HypreJacobi,
+                                        .relative_tol = 1.0e-6,
+                                        .absolute_tol = 1.0e-16,
+                                        .max_iterations = 500,
+                                        .print_level = 0};
+
+  auto equation_solver = std::make_unique<EquationSolver>(nonlin_solver_options, cg_solver_options, pmesh.GetComm());
 
   SolidMechanics<p, dim> solid(std::move(equation_solver), solid_mechanics::default_quasistatic_options, "solid",
-                               mesh_tag);
+                               mesh_tag, std::vector<std::string>{}, 0, 0.0, false, true, 2);
 
   solid_mechanics::NeoHookean mat{.density = 1.0, .K = 1.0, .G = 1.0};