Merge pull request geodynamics#2086 from MFraters/make_spd_function_parameter_available

Make SPD safety factor available as a parameter.

Author: MFraters

include/aspect/newton.h

@@ -186,6 +186,11 @@ namespace aspect
        */
       double get_maximum_linear_stokes_solver_tolerance();
 
+      /**
+       * Get the SPD safety factor.
+       */
+      double get_SPD_safety_factor() const;
+
       /**
        * get a std::string describing the stabilization type used for the preconditioner.
        */
@@ -220,6 +225,7 @@ namespace aspect
       unsigned int        max_newton_line_search_iterations;
       bool                use_newton_residual_scaling_method;
       double              maximum_linear_stokes_solver_tolerance;
+      double              SPD_safety_factor;
   };
 
   namespace Assemblers

include/aspect/utilities.h

@@ -989,7 +989,7 @@ namespace aspect
     double compute_spd_factor(const double eta,
                               const SymmetricTensor<2,dim> &strain_rate,
                               const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
-                              const double safety_factor);
+                              const double SPD_safety_factor);
 
     /**
      * Converts an array of size dim to a Point of size dim.

source/simulator/assemblers/newton_stokes.cc

@@ -120,10 +120,10 @@ namespace aspect
               const SymmetricTensor<2,dim> strain_rate = scratch.material_model_inputs.strain_rate[q];
 
               typename NewtonHandler<dim>::NewtonStabilization preconditioner_stabilization = this->get_newton_handler().get_preconditioner_stabilization();
-              // todo: make this 0.9 into a global input parameter
+
               // use the spd factor when the stabilization is PD or SPD
               const double alpha = (preconditioner_stabilization | NewtonHandler<dim>::NewtonStabilization::PD) != NewtonHandler<dim>::NewtonStabilization::none ?
-                                   Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9)
+                                   Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, this->get_newton_handler().get_SPD_safety_factor())
                                    :
                                    1;
 
@@ -307,10 +307,10 @@ namespace aspect
                   const SymmetricTensor<2,dim> viscosity_derivative_wrt_strain_rate = derivatives->viscosity_derivative_wrt_strain_rate[q];
                   const double viscosity_derivative_wrt_pressure = derivatives->viscosity_derivative_wrt_pressure[q];
                   typename NewtonHandler<dim>::NewtonStabilization velocity_block_stabilization = this->get_newton_handler().get_velocity_block_stabilization();
-                  // todo: make this 0.9 into a global input parameter
+
                   // use the spd factor when the stabilization is PD or SPD
                   const double alpha =  (velocity_block_stabilization | NewtonHandler<dim>::NewtonStabilization::PD) != NewtonHandler<dim>::NewtonStabilization::none ?
-                                        Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, 0.9)
+                                        Utilities::compute_spd_factor<dim>(eta, strain_rate, viscosity_derivative_wrt_strain_rate, this->get_newton_handler().get_SPD_safety_factor())
                                         :
                                         1;
 

source/simulator/newton.cc

@@ -244,6 +244,14 @@ namespace aspect
     return maximum_linear_stokes_solver_tolerance;
   }
 
+  template <int dim>
+  double
+  NewtonHandler<dim>::
+  get_SPD_safety_factor() const
+  {
+    return SPD_safety_factor;
+  }
+
   template <int dim>
   std::string
   NewtonHandler<dim>::
@@ -331,6 +339,18 @@ namespace aspect
                            "When this parameter is true and the linear solver fails, we try again, but now with SPD stabilization "
                            "for both the preconditioner and the velocity block. The SPD stabilization will remain active untill "
                            "the next timestep, when the default values are restored.");
+
+
+        prm.declare_entry ("SPD safety factor", "0.9",
+                           Patterns::Double (0,1),
+                           "When stabilizing the Newton matrix, we can encounter situations where the coefficient inside the elliptic (top-left) "
+                           "block becomes negative or zero. This coefficient has the form $1+x$ where $x$ can sometimes be smaller than $-1$. In "
+                           "this case, the top-left block of the matrix is no longer positive definite, and both preconditioners and iterative "
+                           "solvers may fail. To prevent this, the stabilization computes an $\alpha$ so that $1+\alpha x$ is never negative. "
+                           "This $\alpha$ is chosen as $1$ if $x\ge -1$, and $\alpha=-\frac 1x$ otherwise. (Note that this always leads to "
+                           "$0\le \alpha \le 1$.)  On the other hand, we also want to stay away from $1+\alpha x=0$, and so modify the choice of "
+                           "$\alpha$ to be $1$ if $x\ge -c$, and $\alpha=-\frac cx$ with a $c$ between zero and one. This way, if $c<1$, we are "
+                           "assured that $1-\alpha x>c$, i.e., bounded away from zero.");
       }
       prm.leave_subsection ();
     }
@@ -375,6 +395,9 @@ namespace aspect
 
         AssertThrow((!DEAL_II_VERSION_GTE(9,0,0) && !use_Newton_failsafe) || DEAL_II_VERSION_GTE(9,0,0),
                     ExcMessage("The failsafe option can't be used with a deal.ii less then 9.0.0."));
+
+        SPD_safety_factor = prm.get_double("SPD safety factor");
+
       }
       prm.leave_subsection ();
     }

source/utilities.cc

@@ -2406,7 +2406,7 @@ namespace aspect
     double compute_spd_factor(const double eta,
                               const SymmetricTensor<2,dim> &strain_rate,
                               const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
-                              const double safety_factor)
+                              const double SPD_safety_factor)
     {
       // if the strain rate is zero, or the derivative is zero, then
       // the exact choice of alpha factor does not matter because the
@@ -2420,16 +2420,12 @@ namespace aspect
       const double one_minus_part = 1 - (contract_a_b / norm_a_b);
       const double denom = one_minus_part * one_minus_part * norm_a_b;
 
-      if (denom == 0)
+      // the case denom == 0 (smallest eigenvalue is zero), should return one,
+      // and it does here, because C_safety * 2.0 * eta is always larger then zero.
+      if (denom <= SPD_safety_factor * 2.0 * eta)
         return 1.0;
       else
-        {
-          const double alpha = (2.0*eta)/denom;
-          if (alpha >= 1.0)
-            return 1.0;
-          else
-            return std::max(0.0,safety_factor*alpha);
-        }
+        return std::max(0.0, SPD_safety_factor * ((2.0 * eta) / denom));
     }