Be safer and more correct when using strain rate

source/material_model/rheology/visco_plastic.cc

@@ -118,10 +118,13 @@ namespace aspect
               stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j];
           }
 
-        // The first time this function is called (first iteration of first time step)
-        // a specified "reference" strain rate is used as the returned value would
-        // otherwise be zero.
-        const bool use_reference_strainrate = (this->get_timestep_number() == 0) &&
+        // 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.
+        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());
 
         double edot_ii;
@@ -142,7 +145,12 @@ namespace aspect
               // Choice of activation volume depends on whether there is an adiabatic temperature
               // gradient used when calculating the viscosity. This allows the same activation volume
               // to be used in incompressible and compressible models.
-              const double temperature_for_viscosity = in.temperature[i] + adiabatic_temperature_gradient_for_viscosity*in.pressure[i];
+              const double temperature_for_viscosity = (this->simulator_is_past_initialization())
+                                                       ?
+                                                       in.temperature[i] + adiabatic_temperature_gradient_for_viscosity*in.pressure[i]
+                                                       :
+                                                       this->get_adiabatic_conditions().temperature(in.position[i]);
+
               AssertThrow(temperature_for_viscosity != 0, ExcMessage(
                             "The temperature used in the calculation of the visco-plastic rheology is zero. "
                             "This is not allowed, because this value is used to divide through. It is probably "

source/material_model/visco_plastic.cc

@@ -220,10 +220,12 @@ namespace aspect
           out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic);
           out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic);
 
-          // Compute the effective viscosity if requested and retrieve whether the material is plastically yielding
+          // Compute the effective viscosity if requested and retrieve whether the material is plastically yielding.
+          // Also always compute the viscosity if additional outputs are requested, because the viscosity is needed
+          // to compute the elastic force term.
           bool plastic_yielding = false;
           IsostrainViscosities isostrain_viscosities;
-          if (in.requests_property(MaterialProperties::viscosity))
+          if (in.requests_property(MaterialProperties::viscosity) || in.requests_property(MaterialProperties::additional_outputs))
             {
               // Currently, the viscosities for each of the compositional fields are calculated assuming
               // isostrain amongst all compositions, allowing calculation of the viscosity ratio.

source/simulator/assembly.cc

@@ -555,6 +555,7 @@ namespace aspect
                                            cell,
                                            this->introspection,
                                            current_linearization_point);
+
     scratch.material_model_inputs.requested_properties
       =
         MaterialModel::MaterialProperties::equation_of_state_properties |

source/simulator/melt.cc

@@ -1506,6 +1506,7 @@ namespace aspect
                                          cell,
                                          this->introspection(),
                                          this->get_current_linearization_point());
+
             // We only need access to the MeltOutputs in p_c_scale().
             material_model_inputs.requested_properties = MaterialModel::MaterialProperties::additional_outputs;
 

tests/continental_extension/screen-output

@@ -31,11 +31,11 @@ Number of mesh deformation degrees of freedom: 462
    Solving upper_crust system ... 9 iterations.
    Solving lower_crust system ... 9 iterations.
    Solving lithospheric_mantle system ... 9 iterations.
-   Initial Newton Stokes residual = 8.14683e+16, v = 8.14683e+16, p = 1.18823e+13
+   Initial Newton Stokes residual = 1.832e+14, v = 1.82814e+14, p = 1.18823e+13
 
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+18 iterations.
-      Relative nonlinear residual (total Newton system) after nonlinear iteration 1: 0.000175212, norm of the rhs: 1.42743e+13
+      Relative nonlinear residual (total Newton system) after nonlinear iteration 1: 0.0779163, norm of the rhs: 1.42743e+13
 
 
    Postprocessing:

tests/visco_plastic_derivatives_2d.prm

@@ -7,6 +7,7 @@ set Dimension                              = 2
 set End time                               = 0
 set Use years in output instead of seconds = true
 set Nonlinear solver scheme                = no Advection, no Stokes
+set Adiabatic surface temperature          = 273
 
 # Model geometry (100x100 km, 10 km spacing)
 subsection Geometry model

tests/visco_plastic_derivatives_2d/screen-output

@@ -5,94 +5,94 @@ Loading shared library <./libvisco_plastic_derivatives_2d.debug.so>
 
 Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter harmonic
 pressure: point = 0, Finite difference = 0, Analytical derivative = 0
-pressure: point = 1, Finite difference = 2.12594e+08, Analytical derivative = 2.12594e+08
-pressure: point = 2, Finite difference = 2.68615e+08, Analytical derivative = 2.68615e+08
+pressure: point = 1, Finite difference = 2.45538e+09, Analytical derivative = 2.45538e+09
+pressure: point = 2, Finite difference = 1.01489e+09, Analytical derivative = 1.01489e+09
 pressure: point = 3, Finite difference = 0, Analytical derivative = 0
 pressure: point = 4, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 0, component = 0, Finite difference = 1.14846e+34, Analytical derivative = 1.14846e+34
-strain-rate: point = 1, component = 0, Finite difference = 6.58655e+26, Analytical derivative = 6.58655e+26
-strain-rate: point = 2, component = 0, Finite difference = 1.23761e+27, Analytical derivative = 1.23761e+27
+strain-rate: point = 0, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 0, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 0, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 0, Finite difference = 8.83884e+32, Analytical derivative = 8.83884e+32
-strain-rate: point = 0, component = 1, Finite difference = -1.14846e+34, Analytical derivative = -1.14846e+34
-strain-rate: point = 1, component = 1, Finite difference = -6.58657e+26, Analytical derivative = -6.58657e+26
-strain-rate: point = 2, component = 1, Finite difference = -1.23761e+27, Analytical derivative = -1.23761e+27
+strain-rate: point = 4, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 1, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 1, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 1, Finite difference = -8.83883e+32, Analytical derivative = -8.83883e+32
-strain-rate: point = 0, component = 2, Finite difference = -1.27606e+33, Analytical derivative = -1.27606e+33
-strain-rate: point = 1, component = 2, Finite difference = 1.79805e+28, Analytical derivative = 1.79805e+28
-strain-rate: point = 2, component = 2, Finite difference = -1.1251e+28, Analytical derivative = -1.1251e+28
+strain-rate: point = 4, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 2, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 2, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 2, Finite difference = -4.41942e+32, Analytical derivative = -4.41942e+32
+strain-rate: point = 4, component = 2, Finite difference = 0, Analytical derivative = 0
 OK
 
 Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter geometric
 pressure: point = 0, Finite difference = 0, Analytical derivative = 0
-pressure: point = 1, Finite difference = 2.12594e+08, Analytical derivative = 2.12594e+08
-pressure: point = 2, Finite difference = 2.68615e+08, Analytical derivative = 2.68615e+08
+pressure: point = 1, Finite difference = 2.45538e+09, Analytical derivative = 2.45538e+09
+pressure: point = 2, Finite difference = 1.01489e+09, Analytical derivative = 1.01489e+09
 pressure: point = 3, Finite difference = 0, Analytical derivative = 0
 pressure: point = 4, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 0, component = 0, Finite difference = 1.14845e+34, Analytical derivative = 1.14846e+34
-strain-rate: point = 1, component = 0, Finite difference = 6.58587e+26, Analytical derivative = 6.58655e+26
-strain-rate: point = 2, component = 0, Finite difference = 1.23759e+27, Analytical derivative = 1.23761e+27
+strain-rate: point = 0, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 0, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 0, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 0, Finite difference = 8.83884e+32, Analytical derivative = 8.83884e+32
-strain-rate: point = 0, component = 1, Finite difference = -1.14846e+34, Analytical derivative = -1.14846e+34
-strain-rate: point = 1, component = 1, Finite difference = -6.58596e+26, Analytical derivative = -6.58657e+26
-strain-rate: point = 2, component = 1, Finite difference = -1.23785e+27, Analytical derivative = -1.23761e+27
+strain-rate: point = 4, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 1, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 1, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 1, Finite difference = -8.83883e+32, Analytical derivative = -8.83883e+32
-strain-rate: point = 0, component = 2, Finite difference = -1.27607e+33, Analytical derivative = -1.27606e+33
-strain-rate: point = 1, component = 2, Finite difference = 1.79805e+28, Analytical derivative = 1.79805e+28
-strain-rate: point = 2, component = 2, Finite difference = -1.1251e+28, Analytical derivative = -1.1251e+28
+strain-rate: point = 4, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 2, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 2, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 2, Finite difference = -4.41942e+32, Analytical derivative = -4.41942e+32
+strain-rate: point = 4, component = 2, Finite difference = 0, Analytical derivative = 0
 OK
 
 Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter arithmetic
 pressure: point = 0, Finite difference = 0, Analytical derivative = 0
-pressure: point = 1, Finite difference = 2.12594e+08, Analytical derivative = 2.12594e+08
-pressure: point = 2, Finite difference = 2.68615e+08, Analytical derivative = 2.68615e+08
+pressure: point = 1, Finite difference = 2.45538e+09, Analytical derivative = 2.45538e+09
+pressure: point = 2, Finite difference = 1.01489e+09, Analytical derivative = 1.01489e+09
 pressure: point = 3, Finite difference = 0, Analytical derivative = 0
 pressure: point = 4, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 0, component = 0, Finite difference = 1.14846e+34, Analytical derivative = 1.14846e+34
-strain-rate: point = 1, component = 0, Finite difference = 6.58655e+26, Analytical derivative = 6.58655e+26
-strain-rate: point = 2, component = 0, Finite difference = 1.23761e+27, Analytical derivative = 1.23761e+27
+strain-rate: point = 0, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 0, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 0, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 0, Finite difference = 8.83884e+32, Analytical derivative = 8.83884e+32
-strain-rate: point = 0, component = 1, Finite difference = -1.14846e+34, Analytical derivative = -1.14846e+34
-strain-rate: point = 1, component = 1, Finite difference = -6.58657e+26, Analytical derivative = -6.58657e+26
-strain-rate: point = 2, component = 1, Finite difference = -1.23761e+27, Analytical derivative = -1.23761e+27
+strain-rate: point = 4, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 1, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 1, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 1, Finite difference = -8.83883e+32, Analytical derivative = -8.83883e+32
-strain-rate: point = 0, component = 2, Finite difference = -1.27606e+33, Analytical derivative = -1.27606e+33
-strain-rate: point = 1, component = 2, Finite difference = 1.79805e+28, Analytical derivative = 1.79805e+28
-strain-rate: point = 2, component = 2, Finite difference = -1.1251e+28, Analytical derivative = -1.1251e+28
+strain-rate: point = 4, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 2, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 2, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 2, Finite difference = -4.41942e+32, Analytical derivative = -4.41942e+32
+strain-rate: point = 4, component = 2, Finite difference = 0, Analytical derivative = 0
 OK
 
 Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter maximum composition
 pressure: point = 0, Finite difference = 0, Analytical derivative = 0
-pressure: point = 1, Finite difference = 2.12594e+08, Analytical derivative = 2.12594e+08
-pressure: point = 2, Finite difference = 2.68615e+08, Analytical derivative = 2.68615e+08
+pressure: point = 1, Finite difference = 2.45538e+09, Analytical derivative = 2.45538e+09
+pressure: point = 2, Finite difference = 1.01489e+09, Analytical derivative = 1.01489e+09
 pressure: point = 3, Finite difference = 0, Analytical derivative = 0
 pressure: point = 4, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 0, component = 0, Finite difference = 1.14846e+34, Analytical derivative = 1.14846e+34
-strain-rate: point = 1, component = 0, Finite difference = 6.58655e+26, Analytical derivative = 6.58655e+26
-strain-rate: point = 2, component = 0, Finite difference = 1.23761e+27, Analytical derivative = 1.23761e+27
+strain-rate: point = 0, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 0, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 0, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 0, Finite difference = 8.83884e+32, Analytical derivative = 8.83884e+32
-strain-rate: point = 0, component = 1, Finite difference = -1.14846e+34, Analytical derivative = -1.14846e+34
-strain-rate: point = 1, component = 1, Finite difference = -6.58657e+26, Analytical derivative = -6.58657e+26
-strain-rate: point = 2, component = 1, Finite difference = -1.23761e+27, Analytical derivative = -1.23761e+27
+strain-rate: point = 4, component = 0, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 1, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 1, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 1, Finite difference = -8.83883e+32, Analytical derivative = -8.83883e+32
-strain-rate: point = 0, component = 2, Finite difference = -1.27606e+33, Analytical derivative = -1.27606e+33
-strain-rate: point = 1, component = 2, Finite difference = 1.79805e+28, Analytical derivative = 1.79805e+28
-strain-rate: point = 2, component = 2, Finite difference = -1.1251e+28, Analytical derivative = -1.1251e+28
+strain-rate: point = 4, component = 1, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 0, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 1, component = 2, Finite difference = 0, Analytical derivative = 0
+strain-rate: point = 2, component = 2, Finite difference = 0, Analytical derivative = 0
 strain-rate: point = 3, component = 2, Finite difference = 0, Analytical derivative = 0
-strain-rate: point = 4, component = 2, Finite difference = -4.41942e+32, Analytical derivative = -4.41942e+32
+strain-rate: point = 4, component = 2, Finite difference = 0, Analytical derivative = 0
 OK
 Number of active cells: 100 (on 1 levels)
 Number of degrees of freedom: 2,767 (882+121+441+441+441+441)