From ebc873f05e8a42eb5f7f0409560dfc8a9333438e Mon Sep 17 00:00:00 2001
From: Yimin Jin <jinym@cea-igp.ac.cn>
Date: Mon, 3 Jun 2024 18:05:15 +0800
Subject: [PATCH] generalize the SPD factor

---
 include/aspect/utilities.h | 10 +++-------
 source/utilities.cc        | 23 ++++++++++-------------
 2 files changed, 13 insertions(+), 20 deletions(-)

diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h
index c622377d28a..79350df813f 100644
--- a/include/aspect/utilities.h
+++ b/include/aspect/utilities.h
@@ -841,17 +841,13 @@ namespace aspect
      * this in the Newton paper (Fraters et al., Geophysical Journal
      * International, 2019) where we derived a factor of
      * $\frac{2\eta(\varepsilon(\mathbf u))}{\left[1-\frac{b:a}{\|a\| \|b\|} \right]^2\|a\|\|b\|}$,
-     * which we reset to a maximum of one, and if it is smaller then one,
+     * which we reset to a maximum of one, and if it is smaller than one,
      * a safety_factor scales the value to make sure that 1-alpha won't get to
      * close to zero. However, as later pointed out by Yimin Jin, the computation
      * is wrong, see https://github.com/geodynamics/aspect/issues/5555. Instead,
      * the function now computes the factor as
-     * $(2 \eta) / (a:b + b:a)$, again capped at a maximal value of 1,
-     * and using a safety factor from below.
-     *
-     * In practice, $a$ and $b$ are almost always parallel to each other,
-     * and $a:b + b:a = 2a:b$, in which case one can drop the factor
-     * of $2$ everywhere in the computations.
+     * $\frac{2\eta(\varepsilon(\mathbf u))}{\left[1-\frac{b:a}{\|a\| \|b\|} \right]\|a\|\|b\|}$,
+     * again capped at a maximal value of 1, and using a safety factor from below.
      */
     template <int dim>
     double compute_spd_factor(const double eta,
diff --git a/source/utilities.cc b/source/utilities.cc
index a3c65bc4331..2f9ce449010 100644
--- a/source/utilities.cc
+++ b/source/utilities.cc
@@ -2546,28 +2546,25 @@ namespace aspect
 
       // The factor in the Newton matrix is going to be of the form
       //   2*eta I + (a \otimes b + b \otimes a)
-      // where a=strain_rate and b=dviscosities_dstrain_rate.
-      //
-      // If a,b are parallel, this simplifies to
-      //   [2*eta + 2 a:b] I =  2 [eta + a:b] I
+      // where a=strain_rate and b=dviscosities_dstrain_rate,
       // and we need to make sure that
-      //   [eta + alpha a:b] > (1-safety_factor)*eta
+      //   eta + alpha (a:b - |a||b|) / 2 >= (1-safety_factor)*eta
       // by choosing alpha appropriately.
 
       // So, first check: If
-      //   [eta + a:b] > (1-safety_factor)*eta
+      //   alpha (|a||b| - a:b) <= 2*safety_factor*eta
       // is already satisfied, then we can choose alpha=1
-      const double a_colon_b = strain_rate * dviscosities_dstrain_rate;
-      if (eta + a_colon_b > eta * (1. - SPD_safety_factor))
+      const double norm_a_b = std::sqrt((strain_rate * strain_rate) * (dviscosities_dstrain_rate * dviscosities_dstrain_rate));
+      const double contract_a_b = strain_rate * dviscosities_dstrain_rate;
+      const double denominator = norm_a_b - contract_a_b;
+
+      if (denominator <= 2.0 * SPD_safety_factor * eta)
         return 1.0;
       else
         {
           // Otherwise solve the equation above for alpha, which yields
-          //   a:b = -safety_factor*eta / a:b
-          // This can only ever happen if a:b < 0, so we get
-          //   a:b = safety_factor * abs(eta / a:b)
-          Assert (a_colon_b < 0, ExcInternalError());
-          return SPD_safety_factor * std::abs(eta / a_colon_b);
+          //   alpha = 2*safety_factor*eta / (|a||b| - a:b)
+          return 2.0 * SPD_safety_factor * eta / denominator;
         }
     }