Be safer and more correct when using strain rate

geodynamics · Jul 9, 2023 · 1b37bf5 · 1b37bf5
1 parent dbd7c1c
commit 1b37bf5
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 6 deletions.
diff --git a/source/material_model/rheology/visco_plastic.cc b/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;

diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
@@ -218,9 +218,11 @@ namespace aspect
           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
+          // 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.

diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
@@ -531,7 +531,7 @@ namespace aspect
                                            cell,
                                            this->introspection,
                                            current_linearization_point,
-                                           need_viscosity);
+                                           /*compute_strain_rate=*/ true);
     scratch.material_model_inputs.requested_properties
       =
         MaterialModel::MaterialProperties::equation_of_state_properties |