Add backtracking to ElementJacobianDamper

modules/solid_mechanics/doc/content/source/dampers/ElementJacobianDamper.md

@@ -4,6 +4,10 @@
 
 This damper limits the change in the Jacobians of elements. The damper becomes active if the relative change in Jacobian of any element exceeds the user defined maximum. This damper must be run on the displaced mesh.
 
+Set [!param](/Dampers/ElementJacobianDamper/use_backtracking) to true to have the damper
+iteratively cut back a trial nonlinear update when probing. Without backtracking, a trial update
+that creates a degenerate element map will immediately raise an exception.
+
 !syntax parameters /Dampers/ElementJacobianDamper
 
 !syntax inputs /Dampers/ElementJacobianDamper

modules/solid_mechanics/include/dampers/ElementJacobianDamper.h

@@ -36,6 +36,15 @@ class ElementJacobianDamper : public GeneralDamper
   virtual Real computeDamping(const NumericVector<Number> & /* solution */,
                               const NumericVector<Number> & update) override;
 
+  /**
+   * Probe a damped update on the displaced mesh and report the maximum relative
+   * change in JxW. Returns false if the trial update produces a degenerate map.
+   */
+  bool probeDamping(const NumericVector<Number> & update,
+                    Real damping,
+                    Real & max_difference,
+                    std::string & error_message);
+
 protected:
   /// Thread ID
   THREAD_ID _tid;
@@ -67,4 +76,13 @@ class ElementJacobianDamper : public GeneralDamper
 
   /// Maximum allowed relative increment in Jacobian
   const Real _max_jacobian_diff;
+
+  /// Whether to backtrack the trial update when probing would otherwise create a degenerate map
+  const bool _use_backtracking;
+
+  /// Multiplicative cutback applied during backtracking
+  const Real _backtrack_factor;
+
+  /// Maximum number of backtracking attempts
+  const unsigned int _max_backtrack_steps;
 };

modules/solid_mechanics/src/dampers/ElementJacobianDamper.C

@@ -16,10 +16,24 @@
 #include "libmesh/quadrature.h" // _qrule
 
 // C++
+#include <algorithm>
 #include <cstring> // for "Jacobian" exception test
+#include <limits>
 
 registerMooseObject("SolidMechanicsApp", ElementJacobianDamper);
 
+namespace
+{
+void
+restoreNodes(Elem & elem, const std::vector<Point> & point_copies)
+{
+  mooseAssert(elem.n_nodes() == point_copies.size(), "Node restore cache is the wrong size");
+
+  for (const auto i : make_range(elem.n_nodes()))
+    elem.node_ref(i) = point_copies[i];
+}
+}
+
 InputParameters
 ElementJacobianDamper::validParams()
 {
@@ -31,6 +45,20 @@ ElementJacobianDamper::validParams()
       "max_increment",
       0.1,
       "The maximum permissible relative increment in the Jacobian per Newton iteration");
+  params.addParam<bool>(
+      "use_backtracking",
+      false,
+      "If true, iteratively cut back the probed Newton update when a full trial update would "
+      "create a degenerate displaced element map. The accepted trial must also satisfy "
+      "max_increment.");
+  params.addRangeCheckedParam<Real>(
+      "backtrack_factor",
+      0.5,
+      "backtrack_factor > 0 & backtrack_factor < 1",
+      "Multiplicative cutback applied to the trial damping during ElementJacobianDamper "
+      "backtracking.");
+  params.addParam<unsigned int>(
+      "max_backtrack_steps", 12, "Maximum number of ElementJacobianDamper backtracking attempts.");
   params.set<bool>("use_displaced_mesh") = true;
   return params;
 }
@@ -43,7 +71,10 @@ ElementJacobianDamper::ElementJacobianDamper(const InputParameters & parameters)
     _JxW(_assembly.JxW()),
     _fe_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
     _displaced_problem(_fe_problem.getDisplacedProblem()),
-    _max_jacobian_diff(parameters.get<Real>("max_increment"))
+    _max_jacobian_diff(getParam<Real>("max_increment")),
+    _use_backtracking(getParam<bool>("use_backtracking")),
+    _backtrack_factor(getParam<Real>("backtrack_factor")),
+    _max_backtrack_steps(getParam<unsigned int>("max_backtrack_steps"))
 {
   if (_displaced_problem == NULL)
     mooseError("ElementJacobianDamper: Must use displaced problem");
@@ -69,85 +100,184 @@ Real
 ElementJacobianDamper::computeDamping(const NumericVector<Number> & /* solution */,
                                       const NumericVector<Number> & update)
 {
-  // Maximum difference in the Jacobian for this Newton iteration
-  Real max_difference = 0.0;
+  auto probe_trial =
+      [&](const Real trial_damping, Real & max_difference, std::string & error_message)
+  {
+    const bool trial_nondegenerate_local =
+        probeDamping(update, trial_damping, max_difference, error_message);
+
+    _communicator.max(max_difference);
+
+    const unsigned int invalid_local = trial_nondegenerate_local ? 0 : 1;
+    unsigned int invalid = invalid_local;
+    _communicator.max(invalid);
+
+    if (invalid)
+    {
+      processor_id_type first_invalid_processor =
+          invalid_local ? _communicator.rank() : std::numeric_limits<processor_id_type>::max();
+      _communicator.min(first_invalid_processor);
+      _communicator.broadcast(error_message, first_invalid_processor);
+    }
+
+    return !invalid;
+  };
+
+  if (!_use_backtracking)
+  {
+    Real max_difference = 0.0;
+    std::string error_message;
+
+    PARALLEL_TRY
+    {
+      if (!probe_trial(/*trial_damping=*/1.0, max_difference, error_message))
+        _fe_problem.setException(error_message);
+    }
+    PARALLEL_CATCH;
+
+    if (max_difference > _max_jacobian_diff)
+      return _max_jacobian_diff / max_difference;
+
+    return 1.0;
+  }
+
+  const Real minimum_trial_damping = std::max(_min_damping, std::numeric_limits<Real>::epsilon());
+  Real damping = 1.0;
+  std::string error_message;
+
+  // Try 1 step, then at most max_backtrack_steps extra steps
+  for (unsigned int step = 0; step <= _max_backtrack_steps; ++step)
+  {
+    Real max_difference = 0.0;
+    bool trial_nondegenerate = false;
+
+    PARALLEL_TRY { trial_nondegenerate = probe_trial(damping, max_difference, error_message); }
+    PARALLEL_CATCH;
+
+    if (trial_nondegenerate && max_difference <= _max_jacobian_diff)
+      return damping;
+
+    if (damping <= minimum_trial_damping || step == _max_backtrack_steps)
+    {
+      if (!trial_nondegenerate)
+        _fe_problem.setException(error_message.empty()
+                                     ? "ElementJacobianDamper could not find a nondegenerate "
+                                       "trial update."
+                                     : error_message);
+      else
+        _fe_problem.setException("ElementJacobianDamper could not reduce the relative Jacobian "
+                                 "increment below max_increment without driving the damping "
+                                 "factor below a usable threshold.");
+
+      // The exception flag will abort the nonlinear step after the damper returns, so return the
+      // last trial damping only to satisfy the interface.
+      return damping;
+    }
+
+    if (!trial_nondegenerate)
+      damping *= _backtrack_factor;
+    else
+      damping =
+          std::min(damping * _backtrack_factor, damping * _max_jacobian_diff / max_difference);
+  }
+
+  return damping;
+}
+
+bool
+ElementJacobianDamper::probeDamping(const NumericVector<Number> & update,
+                                    const Real damping,
+                                    Real & max_difference,
+                                    std::string & error_message)
+{
   MooseArray<Real> JxW_displaced;
 
   // Vector for storing the original node coordinates
   std::vector<Point> point_copies;
 
-  PARALLEL_TRY
+  max_difference = 0.0;
+  error_message.clear();
+
+  try
   {
-    try
+    // Loop over elements in the mesh
+    for (auto & current_elem : _mesh->getMesh().active_local_element_ptr_range())
     {
-      // Loop over elements in the mesh
-      for (auto & current_elem : _mesh->getMesh().active_local_element_ptr_range())
-      {
-        point_copies.clear();
-        point_copies.reserve(current_elem->n_nodes());
+      point_copies.clear();
+      point_copies.reserve(current_elem->n_nodes());
 
-        // Displace nodes with current Newton increment
-        for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
-        {
-          Node & displaced_node = current_elem->node_ref(i);
+      // Displace nodes with the trial Newton increment
+      for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
+      {
+        Node & displaced_node = current_elem->node_ref(i);
 
-          point_copies.push_back(displaced_node);
+        point_copies.push_back(displaced_node);
 
-          for (unsigned int j = 0; j < _ndisp; ++j)
-          {
-            dof_id_type disp_dof_num =
-                displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
-            displaced_node(j) += update(disp_dof_num);
-          }
+        for (unsigned int j = 0; j < _ndisp; ++j)
+        {
+          const dof_id_type disp_dof_num =
+              displaced_node.dof_number(_sys.number(), _disp_var[j]->number(), 0);
+          displaced_node(j) += damping * update(disp_dof_num);
         }
+      }
 
-        // Reinit element to compute Jacobian of displaced element
+      try
+      {
+        // Reinit element to compute Jacobian of the trial displaced element
         _assembly.reinit(current_elem);
         JxW_displaced = _JxW;
+      }
+      catch (const std::exception & e)
+      {
+        restoreNodes(*current_elem, point_copies);
 
-        // Un-displace nodes
-        for (unsigned int i = 0; i < current_elem->n_nodes(); ++i)
+        // Degenerate-map failures mean this trial damping is too aggressive
+        if (strstr(e.what(), "Jacobian") || strstr(e.what(), "singular") ||
+            strstr(e.what(), "det != 0"))
         {
-          Node & displaced_node = current_elem->node_ref(i);
-
-          for (unsigned int j = 0; j < _ndisp; ++j)
-            displaced_node(j) = point_copies[i](j);
+          error_message = e.what();
+          return false;
         }
 
-        // Reinit element to compute Jacobian before displacement
-        _assembly.reinit(current_elem);
+        throw;
+      }
+
+      restoreNodes(*current_elem, point_copies);
 
-        for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
+      // Reinit element to compute Jacobian before displacement
+      try
+      {
+        _assembly.reinit(current_elem);
+      }
+      catch (const std::exception & e)
+      {
+        if (strstr(e.what(), "Jacobian") || strstr(e.what(), "singular") ||
+            strstr(e.what(), "det != 0"))
         {
-          Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / _JxW[qp];
-          if (diff > max_difference)
-            max_difference = diff;
+          error_message = e.what();
+          return false;
         }
 
-        JxW_displaced.release();
-      }
-    }
-    catch (MooseException & e)
-    {
-      _fe_problem.setException(e.what());
-    }
-    catch (std::exception & e)
-    {
-      // Continue if we find a libMesh degenerate map exception, but
-      // just throw for any real error
-      if (!strstr(e.what(), "Jacobian") && !strstr(e.what(), "singular") &&
-          !strstr(e.what(), "det != 0"))
         throw;
+      }
+
+      for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
+      {
+        const Real denominator = std::max(std::abs(_JxW[qp]), libMesh::TOLERANCE);
+        const Real diff = std::abs(JxW_displaced[qp] - _JxW[qp]) / denominator;
+        if (diff > max_difference)
+          max_difference = diff;
+      }
 
-      _fe_problem.setException(e.what());
+      JxW_displaced.release();
     }
   }
-  PARALLEL_CATCH;
-
-  _communicator.max(max_difference);
-
-  if (max_difference > _max_jacobian_diff)
-    return _max_jacobian_diff / max_difference;
+  catch (const MooseException & e)
+  {
+    error_message = e.what();
+    _fe_problem.setException(e.what());
+    return false;
+  }
 
-  return 1.0;
+  return true;
 }

modules/solid_mechanics/test/tests/jacobian_damper/tests

@@ -38,4 +38,12 @@
     expect_exit_code = 0
     expect_err = 'negative Jacobian'
   []
+  [backtrack_negative_jacobian]
+    type = RunApp
+    input = 'exception_handling.i'
+    cli_args = 'Dampers/jac/use_backtracking=true'
+    max_parallel = 1
+    expect_out = 'Damping factor:'
+    requirement = "The element Jacobian damper shall optionally backtrack a trial nonlinear update to avoid generating a degenerate element map."
+  []
 []