Merge pull request #5981 from anne-glerum/use_min_strain_rate

Only use reference strain rate during the first nonlinear iteration in viscoplastic material model

source/material_model/rheology/visco_plastic.cc

@@ -118,14 +118,14 @@ namespace aspect
               stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j];
           }
 
-        // Use a specified "reference" strain rate if the strain rate is not yet available,
-        // or close to zero. This is to avoid division by zero.
+        // Use a specified "reference" strain rate if the strain rate is not yet available
+        // during the very first nonlinear iteration or before. This is to avoid division by zero and
+        // to have a better estimate of the resulting viscosity during this time.
+        // During later iterations and timesteps the strain rate is capped by a minimum value min_strain_rate.
         const bool use_reference_strainrate = this->simulator_is_past_initialization() == false
                                               ||
                                               (this->get_timestep_number() == 0 &&
-                                               this->get_nonlinear_iteration() == 0)
-                                              ||
-                                              (in.strain_rate[i].norm() <= std::numeric_limits<double>::min());
+                                               this->get_nonlinear_iteration() == 0);
 
         double edot_ii;
         if (use_reference_strainrate)
@@ -262,9 +262,7 @@ namespace aspect
               {
                 const std::vector<double> &elastic_shear_moduli = elastic_rheology.get_elastic_shear_moduli();
 
-                if (use_reference_strainrate == true)
-                  effective_edot_ii = ref_strain_rate;
-                else
+                if (use_reference_strainrate == false)
                   {
                     // Overwrite effective_edot_ii with a value that includes a term that accounts for
                     // elastic stress arising from a previous time step.
@@ -650,6 +648,8 @@ namespace aspect
         // Reference and minimum/maximum values
         min_strain_rate = prm.get_double("Minimum strain rate");
         ref_strain_rate = prm.get_double("Reference strain rate");
+        AssertThrow(ref_strain_rate >= min_strain_rate,
+                    ExcMessage("The reference strain rate for the viscoplastic material model should be larger than the minimum strain rate."));
 
         // Retrieve the list of composition names
         std::vector<std::string> compositional_field_names = this->introspection().get_composition_names();

tests/continental_extension/screen-output

@@ -34,11 +34,11 @@ Continuing to the next timestep even though solution is not fully converged.
    Solving upper_crust system ... 9 iterations.
    Solving lower_crust system ... 9 iterations.
    Solving lithospheric_mantle system ... 9 iterations.
-   Initial Newton Stokes residual = 1.57637e+14, v = 1.57164e+14, p = 1.21967e+13
+   Initial Newton Stokes residual = 7.74397e+16, v = 7.74397e+16, p = 1.21967e+13
 
    Solving Stokes system... 0+21 iterations.
       Newton system information: Norm of the rhs: 2.35285e+13
-      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.149257
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.00030383
 
 
 WARNING: The nonlinear solver in the current timestep failed to converge.

tests/visco_plastic_yield_strain_weakening_defect_corr_stokes/screen-output

@@ -26,11 +26,11 @@ Continuing to the next timestep even though solution is not fully converged.
 *** Timestep 1:  t=1 years, dt=1 years
    Solving temperature system... 0 iterations.
    Solving total_strain system ... 10 iterations.
-   Initial Newton Stokes residual = 1.8223e+13, v = 1.81808e+13, p = 1.24042e+12
+   Initial Newton Stokes residual = 1.85485e+13, v = 1.85069e+13, p = 1.24042e+12
 
    Solving Stokes system... 0+6 iterations.
       Newton system information: Norm of the rhs: 1.80069e+09
-      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 9.8814e-05
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 9.70804e-05
 
 
 WARNING: The nonlinear solver in the current timestep failed to converge.