make deviator() and second_invariant() consistent with plane strain assumption in 2D

Author: YiminJin

include/aspect/utilities.h

@@ -1407,6 +1407,29 @@ namespace aspect
       // Declare the existence of a specialization:
       template <>
       const Tensor<3,3> &levi_civita<3>();
+
+      /**
+       * Compute the deviator of a symmetric tensor. Unlike the function with the same name in
+       * deal.II, this function is consistent with the plane strain assumption when dim = 2.
+       * Specifically, the deviator of the stress tensor $\mathbf\tau$ in 2D is given by 
+       * $\text{dev}(\mathbf\tau) = \mathbf\tau - \frac{1}{3}\text{trace}(\mathbf\tau)\mathbf 1$
+       * under the plane strain assumption.
+       */
+      template <int dim>
+      SymmetricTensor<2,dim> 
+      deviator(const SymmetricTensor<2,dim> &input);
+
+      /**
+       * Compute the second invariant of a deviatoric tensor of rank 2. Unlike the function
+       * named "second_invariant" in deal.II, this function is consistent with the plane
+       * strain assumption when dim = 2. Specifically, the second invariant of the deviatoric
+       * stress tensor $\tau_{II}$ in 2D is given by $\tau_{II} = -\frac{1}{2}(\tau_{11}^2 +
+       * \tau_{22}^2 + \tau_{33}^2 + 2\tau_{12}^2) = -\frac{1}{2}[\tau_{11}^2 + \tau_{22}^2 +
+       * (\tau_{11} + \tau_{22})^2 + 2\tau_{12}^2]$ under the plane strain assumption.
+       */
+      template <int dim>
+      double
+      deviatoric_tensor_inv2(const SymmetricTensor<2,dim> &input);
     }
 
   }

source/heating_model/shear_heating.cc

@@ -76,7 +76,7 @@ namespace aspect
                                                                                          drucker_prager_parameters);
 
               // Scale the stress accordingly.
-              const double deviatoric_stress = 2 * material_model_outputs.viscosities[q] * std::sqrt(std::fabs(second_invariant(deviatoric_strain_rate)));
+              const double deviatoric_stress = 2 * material_model_outputs.viscosities[q] * std::sqrt(std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(deviatoric_strain_rate)));
               double scaling_factor = 1.0;
               if (deviatoric_stress > 0.0)
                 scaling_factor = std::min(yield_stress / deviatoric_stress, 1.0);

source/material_model/drucker_prager.cc

@@ -53,7 +53,7 @@ namespace aspect
               Assert(std::isfinite(in.strain_rate[i].norm()),
                      ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
                                 "not filled by the caller."));
-              const SymmetricTensor<2,dim> strain_rate_deviator = deviator(in.strain_rate[i]);
+              const SymmetricTensor<2,dim> strain_rate_deviator = Utilities::Tensors::deviator(in.strain_rate[i]);
 
               // For the very first time this function is called
               // (the first iteration of the first timestep), this function is called
@@ -64,8 +64,8 @@ namespace aspect
               // In later iterations and timesteps we calculate the second moment
               // invariant of the deviatoric strain rate tensor.
               // This is equal to the negative of the second principle
-              // invariant of the deviatoric strain rate (calculated with the function second_invariant),
-              // as shown in Appendix A of Zienkiewicz and Taylor (Solid Mechanics, 2000).
+              // invariant of the deviatoric strain rate, as shown in Appendix A of 
+              // Zienkiewicz and Taylor (Solid Mechanics, 2000).
               //
               // The negative of the second principle invariant is equal to 0.5 e_dot_dev_ij e_dot_dev_ji,
               // where e_dot_dev is the deviatoric strain rate tensor. The square root of this quantity
@@ -78,7 +78,7 @@ namespace aspect
                                              // initialized. This might mean that we are
                                              // in a unit test, or at least that we can't
                                              // rely on any simulator information
-                                             std::fabs(second_invariant(strain_rate_deviator))
+                                             std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(strain_rate_deviator))
                                              :
                                              // simulator object is available, but we need to treat the
                                              // first time step separately
@@ -88,7 +88,7 @@ namespace aspect
                                               ?
                                               reference_strain_rate * reference_strain_rate
                                               :
-                                              std::fabs(second_invariant(strain_rate_deviator))));
+                                              std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(strain_rate_deviator))));
 
               const double strain_rate_effective = edot_ii_strict;
 

source/material_model/entropy_model.cc

@@ -253,7 +253,7 @@ namespace aspect
                                                                                                       drucker_prager_parameters);
                     }
 
-                  const double strain_rate_effective = std::fabs(second_invariant(deviator(in.strain_rate[i])));
+                  const double strain_rate_effective = std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])));
 
                   if (std::sqrt(strain_rate_effective) >= std::numeric_limits<double>::min())
                     {

source/material_model/grain_size.cc

@@ -199,7 +199,7 @@ namespace aspect
 
       const double strain_rate_dependence = (1.0 - dislocation_creep_exponent[phase_index]) / dislocation_creep_exponent[phase_index];
       const SymmetricTensor<2,dim> shear_strain_rate = strain_rate - 1./dim * trace(strain_rate) * unit_symmetric_tensor<dim>();
-      const double second_strain_rate_invariant = std::sqrt(std::max(-second_invariant(shear_strain_rate), 0.));
+      const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(shear_strain_rate), 0.));
 
       // If the strain rate is zero, the dislocation viscosity is infinity.
       if (second_strain_rate_invariant <= std::numeric_limits<double>::min())
@@ -221,7 +221,7 @@ namespace aspect
         {
           const SymmetricTensor<2,dim> dislocation_strain_rate = diffusion_viscosity
                                                                  / (diffusion_viscosity + dis_viscosity) * shear_strain_rate;
-          const double dislocation_strain_rate_invariant = std::sqrt(std::max(-second_invariant(dislocation_strain_rate), 0.));
+          const double dislocation_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(dislocation_strain_rate), 0.));
 
           dis_viscosity_old = dis_viscosity;
           dis_viscosity = dislocation_creep_prefactor[phase_index]
@@ -540,8 +540,8 @@ namespace aspect
               Assert(std::isfinite(in.strain_rate[i].norm()),
                      ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
                                 "not filled by the caller."));
-              const SymmetricTensor<2,dim> shear_strain_rate = deviator(in.strain_rate[i]);
-              const double second_strain_rate_invariant = std::sqrt(std::max(-second_invariant(shear_strain_rate), 0.));
+              const SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::deviator(in.strain_rate[i]);
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(shear_strain_rate), 0.));
 
               const double adiabatic_temperature = this->get_adiabatic_conditions().is_initialized()
                                                    ?

source/material_model/reaction_model/grain_size_evolution.cc

@@ -227,7 +227,7 @@ namespace aspect
 
               // grain size reduction in dislocation creep regime
               const SymmetricTensor<2,dim> shear_strain_rate = in.strain_rate[i] - 1./dim * trace(in.strain_rate[i]) * unit_symmetric_tensor<dim>();
-              const double second_strain_rate_invariant = std::sqrt(std::max(-second_invariant(shear_strain_rate), 0.));
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(shear_strain_rate), 0.));
 
               const double current_diffusion_viscosity   = diffusion_viscosity(in.temperature[i], adiabatic_temperature, pressures[i], grain_size, second_strain_rate_invariant, phase_indices[i]);
               current_dislocation_viscosity = dislocation_viscosity(in.temperature[i], adiabatic_temperature, pressures[i], in.strain_rate[i], phase_indices[i], current_diffusion_viscosity, current_dislocation_viscosity);

source/material_model/rheology/composite_visco_plastic.cc

@@ -161,7 +161,7 @@ namespace aspect
         // to prevent a division-by-zero, and a floating point exception.
         // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
         // strain rate (often simplified as epsilondot_ii)
-        const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
                                         minimum_strain_rate);
         const double log_edot_ii = std::log(edot_ii);
 
@@ -475,7 +475,7 @@ namespace aspect
         // to prevent a division-by-zero, and a floating point exception.
         // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
         // strain rate (often simplified as epsilondot_ii)
-        const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
                                         minimum_strain_rate);
         const double log_edot_ii = std::log(edot_ii);
 

source/material_model/rheology/diffusion_dislocation.cc

@@ -51,7 +51,7 @@ namespace aspect
         // to prevent a division-by-zero, and a floating point exception.
         // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
         // strain rate (often simplified as epsilondot_ii)
-        const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
                                         min_strain_rate);
         const double log_edot_ii = std::log(edot_ii);
 

source/material_model/rheology/elasticity.cc

@@ -401,7 +401,7 @@ namespace aspect
 
             for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
               {
-                const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]);
+                const SymmetricTensor<2, dim> deviatoric_strain_rate = Utilities::Tensors::deviator(in.strain_rate[i]);
 
                 // Get stress from timestep $t$ rotated and advected into the current
                 // timestep $t+\Delta t_c$ from the compositional fields.
@@ -496,7 +496,7 @@ namespace aspect
         for (unsigned int i = 0; i < in.n_evaluation_points(); ++i)
           {
             const double eta = out.viscosities[i];
-            const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]);
+            const SymmetricTensor<2, dim> deviatoric_strain_rate = Utilities::Tensors::deviator(in.strain_rate[i]);
             const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][stress_start_index],
                                                             &in.composition[i][stress_start_index]+n_independent_components));
             const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][stress_start_index+n_independent_components],
@@ -700,7 +700,7 @@ namespace aspect
 
                 // Compute the total stress at time t.
                 const SymmetricTensor<2, dim>
-                stress_t = 2. * effective_creep_viscosity * deviator(in.strain_rate[i])
+                stress_t = 2. * effective_creep_viscosity * Utilities::Tensors::deviator(in.strain_rate[i])
                            + effective_creep_viscosity / elastic_viscosity * stress_0_t
                            + (1. - timestep_ratio) * (1. - effective_creep_viscosity / elastic_viscosity) * stress_old;
 
@@ -858,7 +858,7 @@ namespace aspect
         const double timestep_ratio = calculate_timestep_ratio();
 
         const SymmetricTensor<2, dim>
-        edot_deviator = deviator(strain_rate) + 0.5 * stress_0_advected / elastic_viscosity
+        edot_deviator = Utilities::Tensors::deviator(strain_rate) + 0.5 * stress_0_advected / elastic_viscosity
                         + 0.5 * (1. - timestep_ratio) * (1.  - creep_viscosity/elastic_viscosity) * stress_old / creep_viscosity;
 
         return edot_deviator;

source/material_model/rheology/strain_dependent.cc

@@ -457,7 +457,7 @@ namespace aspect
                                                          &composition[0] + Tensor<2,dim>::n_independent_components));
               const SymmetricTensor<2,dim> L = symmetrize( strain * transpose(strain) );
 
-              const double strain_ii = std::fabs(second_invariant(L));
+              const double strain_ii = std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(L));
               brittle_weakening = calculate_plastic_weakening(strain_ii, j);
               viscous_weakening = calculate_viscous_weakening(strain_ii, j);
               break;
@@ -633,7 +633,7 @@ namespace aspect
                    ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
                               "not filled by the caller."));
 
-            const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(in.strain_rate[i])), 0.)),
+            const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])), 0.)),
                                             min_strain_rate);
             double delta_e_ii = edot_ii*this->get_timestep();