try again.

geodynamics · Mar 1, 2018 · 9dae3f4 · 9dae3f4
1 parent f506c72
commit 9dae3f4
Showing 1 changed file with 3 additions and 3 deletions.
diff --git a/benchmarks/nonlinear_channel_flow/simple_nonlinear.cc b/benchmarks/nonlinear_channel_flow/simple_nonlinear.cc
@@ -200,8 +200,8 @@ namespace aspect
                   // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
                   // strain rate (often simplified as epsilondot_ii)
                   const SymmetricTensor<2,dim> edot = use_deviator_of_strain_rate ? deviator(in.strain_rate[i]) : in.strain_rate[i];
-                  const double edot_ii_strict = std::max(std::sqrt(edot*edot), std::sqrt(2) * min_strain_rate[c]);
-                  const double edot_ii = 2.0 * std::sqrt(0.5) * edot_ii_strict;
+                  const double edot_ii_strict = std::sqrt(0.5*edot*edot);
+                  const double edot_ii = 2.0 * std::max(edot_ii_strict, min_strain_rate[c]);
 
                   const double stress_exponent_inv = 1/stress_exponent[c];
                   composition_viscosities[c] = std::max(std::min(std::pow(viscosity_prefactor[c],-stress_exponent_inv) * std::pow(edot_ii,stress_exponent_inv-1), max_viscosity[c]), min_viscosity[c]);
@@ -210,7 +210,7 @@ namespace aspect
                     {
                       if (edot_ii_strict > min_strain_rate[c] && composition_viscosities[c] < max_viscosity[c] && composition_viscosities[c] > min_viscosity[c])
                         {
-                          composition_viscosities_derivatives[c] = (stress_exponent_inv-1) * composition_viscosities[c] * (1.0/(edot_ii_strict*edot_ii_strict)) * edot;
+                          composition_viscosities_derivatives[c] = (stress_exponent_inv-1) * composition_viscosities[c] * (1.0/(edot*edot)) * edot;
 
                           if (use_deviator_of_strain_rate == true)
                             composition_viscosities_derivatives[c] = composition_viscosities_derivatives[c] * deviator_tensor<dim>();