Improve documentation

benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/run.sh

@@ -28,8 +28,8 @@ for case in 1 2 3; do
   mkdir -p $dirname
 
   sed  \
-    -e "s/set Function expression = if(x<60e3,.*/    set Function expression = if(x<60e3,$U,-$U);0/g" \
-    -e "s/set Reference viscosity .*/    set Reference viscosity = $background_viscosity/g" \
+    -e "s/set Function expression = if(x<60e3,.*/set Function expression = if(x<60e3,$U,-$U);0/g" \
+    -e "s/set Reference viscosity .*/set Reference viscosity = $background_viscosity/g" \
     -e "s/set Output directory .*/set Output directory = results\/$dirname_base/g" \
     -e "s/set Nonlinear solver scheme .*/set Nonlinear solver scheme = single Advection, iterated Stokes/g" \
     -e "s/set List of output variables .*/set List of output variables = material properties,strain rate/g" \
@@ -49,10 +49,10 @@ for case in 1 2 3; do
         mkdir -p $dirname
 
         sed  \
-          -e "s/set Function expression = if(x<60e3,.*/    set Function expression = if(x<60e3,$U,-$U);0/g" \
-          -e "s/set Reference viscosity .*/    set Reference viscosity = $background_viscosity/g" \
+          -e "s/set Function expression = if(x<60e3,.*/set Function expression = if(x<60e3,$U,-$U);0/g" \
+          -e "s/set Reference viscosity .*/set Reference viscosity = $background_viscosity/g" \
           -e "s/set Output directory .*/set Output directory = results\/$dirname_base/g" \
-          -e "s/set Max pre-Newton nonlinear iterations .*/    set Max pre-Newton nonlinear iterations = $picard_iterations/g" \
+          -e "s/set Max pre-Newton nonlinear iterations .*/set Max pre-Newton nonlinear iterations = $picard_iterations/g" \
           -e "s/set Stabilization preconditioner .*/set Stabilization preconditioner = $stabilization/g" \
           -e "s/set Stabilization velocity block .*/set Stabilization velocity block = $stabilization/g" \
           -e "s/set Maximum linear Stokes solver tolerance .*/set Maximum linear Stokes solver tolerance = $maximum_linear_tolerance/g" \

source/simulator/assemblers/newton_stokes.cc

@@ -178,6 +178,10 @@ namespace aspect
               const typename Newton::Parameters::Stabilization
               preconditioner_stabilization = this->get_newton_handler().parameters.preconditioner_stabilization;
 
+              // use the correct strain rate for the Jacobian
+              // when elasticity is enabled use viscoelastic strain rate
+              // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+              // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
               SymmetricTensor<2,dim> effective_strain_rate = scratch.material_model_inputs.strain_rate[q];
               if (elastic_out != nullptr)
                 effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];
@@ -479,6 +483,10 @@ namespace aspect
                   const Newton::Parameters::Stabilization velocity_block_stabilization
                     = this->get_newton_handler().parameters.velocity_block_stabilization;
 
+                  // use the correct strain rate for the Jacobian
+                  // when elasticity is enabled use viscoelastic strain rate
+                  // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+                  // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
                   SymmetricTensor<2,dim> effective_strain_rate = scratch.material_model_inputs.strain_rate[q];
                   if (elastic_out != nullptr)
                     effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];

source/simulator/stokes_matrix_free.cc

@@ -1424,6 +1424,10 @@ namespace aspect
 
                   for (unsigned int q=0; q<n_q_points; ++q)
                     {
+                      // use the correct strain rate for the Jacobian
+                      // when elasticity is enabled use viscoelastic strain rate
+                      // when stabilization is enabled, use the deviatoric strain rate because the SPD factor
+                      // that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
                       SymmetricTensor<2,dim> effective_strain_rate = in.strain_rate[q];
                       if (elastic_out != nullptr)
                         effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];