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

benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/drucker_prager_compositions.cc

@@ -32,6 +32,7 @@
 #include <aspect/material_model/interface.h>
 #include <aspect/simulator_access.h>
 #include <aspect/newton.h>
+#include <aspect/utilities.h>
 
 namespace aspect
 {
@@ -288,7 +289,7 @@ namespace aspect
               std::vector<SymmetricTensor<2,dim>> composition_viscosities_derivatives(volume_fractions.size());
               std::vector<double> composition_dviscosities_dpressure(volume_fractions.size());
 
-              const SymmetricTensor<2,dim> deviator_strain_rate = use_deviator_of_strain_rate ? deviator(in.strain_rate[i]) : in.strain_rate[i];
+              const SymmetricTensor<2,dim> deviator_strain_rate = use_deviator_of_strain_rate ? Utilities::Tensors::consistent_deviator(in.strain_rate[i]) : in.strain_rate[i];
 
               for (unsigned int c=0; c < volume_fractions.size(); ++c)
                 {
@@ -297,7 +298,9 @@ 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 double edot_ii = compute_second_invariant(deviator_strain_rate, min_strain_rate[c]);
+                  const double edot_ii = sqrt(2.0 * std::abs(use_deviator_of_strain_rate ?
+                                                             Utilities::Tensors::consistent_second_invariant_of_deviatoric_tensor(deviator_strain_rate) :
+                                                             second_invariant(deviator_strain_rate)));
 
                   // Find effective viscosities for each of the individual phases
                   // Viscosities should have same number of entries as compositional fields
@@ -316,15 +319,11 @@ namespace aspect
                             {
                               const double drucker_prager_viscosity = compute_viscosity(edot_ii,pressure,c,prefactor[c],false,min_visc[c],max_visc[c]);
                               const double regularization_adjustment = (ref_visc * ref_visc)
-                                                                       / (ref_visc * ref_visc + 2.0 * ref_visc * drucker_prager_viscosity
-                                                                          + drucker_prager_viscosity * drucker_prager_viscosity);
+                                                                       / ((ref_visc + drucker_prager_viscosity) * (ref_visc + drucker_prager_viscosity));
 
                               composition_viscosities_derivatives[c] = -regularization_adjustment *
                                                                        (drucker_prager_viscosity / (edot_ii * edot_ii)) * deviator_strain_rate;
 
-                              if (use_deviator_of_strain_rate == true)
-                                composition_viscosities_derivatives[c]*deviator_tensor<dim>();
-
                               composition_dviscosities_dpressure[c] = regularization_adjustment *
                                                                       ((dim == 3)
                                                                        ?
@@ -349,9 +348,14 @@ namespace aspect
                               const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component);
 
                               // add a small difference to one independent component of the strain-rate tensor
-                              const SymmetricTensor<2,dim> strain_rate_difference_plus = deviator_strain_rate + std::max(edot_ii, min_strain_rate[c]*min_strain_rate[c]) * finite_difference_accuracy
-                                                                                         * Utilities::nth_basis_for_symmetric_tensors<dim>(component);
-                              const double second_invariant_strain_rate_difference_plus = compute_second_invariant(strain_rate_difference_plus, min_strain_rate[c]);
+                              SymmetricTensor<2,dim> strain_rate_difference_plus = in.strain_rate[i] + std::max(edot_ii, min_strain_rate[c]*min_strain_rate[c]) * finite_difference_accuracy
+                                                                                   * Utilities::nth_basis_for_symmetric_tensors<dim>(component);
+                              if (use_deviator_of_strain_rate)
+                                strain_rate_difference_plus = Utilities::Tensors::consistent_deviator(strain_rate_difference_plus);
+                              const double second_invariant_strain_rate_difference_plus =
+                                std::sqrt(2.0 * std::abs(use_deviator_of_strain_rate ?
+                                                         Utilities::Tensors::consistent_second_invariant_of_deviatoric_tensor(strain_rate_difference_plus) :
+                                                         second_invariant(strain_rate_difference_plus)));
                               const double eta_component_plus = compute_viscosity(second_invariant_strain_rate_difference_plus,pressure,c,prefactor[c],true,min_visc[c],max_visc[c]);
 
                               // compute the difference between the viscosity with and without the strain-rate difference.

doc/modules/changes/20250617_yiminjin

@@ -0,0 +1,7 @@
+Changed: The deviator and second invariant of symmetric tensors are modified to
+be consistent with the plane strain assumption in 2D. This fixes the
+inconsistency between some material models and the Stokes assemblers when
+computing the deviatoric strain rate. It also results in better nonlinear solver
+convergence due to some fixes to the Newton solver.
+<br>
+(Yimin Jin, 2025/06/17)

include/aspect/newton.h

@@ -262,6 +262,18 @@ namespace aspect
         void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
     };
 
+    /**
+     * This class assembles the compressible adjustment to the Newton Stokes preconditioner.
+     */
+    template <int dim>
+    class NewtonStokesCompressiblePreconditioner : public NewtonInterface<dim>
+    {
+      public:
+        void
+        execute (internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
+                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
+    };
+
     /**
      * This class assembles the terms for the matrix and right-hand-side of the incompressible
      * Newton Stokes system for the current cell.

include/aspect/simulator/solver/matrix_free_operators.h

@@ -122,6 +122,11 @@ namespace aspect
        */
       bool apply_stabilization_free_surface_faces;
 
+      /**
+       * If true, average the Newton factors in each cell.
+       */
+      bool average_newton_factors;
+
       /**
        * Table which stores viscosity values for each cell.
        *

include/aspect/utilities.h

@@ -1379,6 +1379,38 @@ namespace aspect
       // Declare the existence of a specialization:
       template <>
       const Tensor<3,3> &levi_civita<3>();
+
+      /**
+       * Compute the deviator of a symmetric tensor. This function is equivalent to
+       * dealii::deviator in 3D, while in 2D it is consistent with the plane strain assumption.
+       * 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.
+       *
+       * It should be noted that the consistent deviator of a consistently deviatoric tensor is
+       * not itself in 2D, which implies that it is invalid to apply this function to a symmetric
+       * tensor more than once.
+       */
+      template <int dim>
+      SymmetricTensor<2,dim>
+      consistent_deviator(const SymmetricTensor<2,dim> &input);
+
+      /**
+       * Compute the second invariant of a deviatoric tensor of rank 2. This function is
+       * equivalent to dealii::second_invariant in 3D, while in 2D it is consistent with the
+       * plane strain assumption. 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.
+       *
+       * It should be noted that this function provides the correct result in 2D only when the
+       * input tensor is deviatoric in the sense of plane strain, which cannot be examined
+       * without extra information (the trace of a plane strain deviatoric tensor is not
+       * zero). Thus, extra care must be taken when using this function.
+       */
+      template <int dim>
+      double
+      consistent_second_invariant_of_deviatoric_tensor(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::consistent_second_invariant_of_deviatoric_tensor(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::consistent_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::consistent_second_invariant_of_deviatoric_tensor(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::consistent_second_invariant_of_deviatoric_tensor(strain_rate_deviator))));
 
               const double strain_rate_effective = edot_ii_strict;
 

source/material_model/entropy_model.cc

@@ -416,7 +416,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::consistent_second_invariant_of_deviatoric_tensor(Utilities::Tensors::consistent_deviator(in.strain_rate[i])));
 
                   if (std::sqrt(strain_rate_effective) >= std::numeric_limits<double>::min())
                     {

source/material_model/grain_size.cc

@@ -198,8 +198,8 @@ 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 SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::consistent_deviator(strain_rate);
+      const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::consistent_second_invariant_of_deviatoric_tensor(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::consistent_second_invariant_of_deviatoric_tensor(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::consistent_deviator(in.strain_rate[i]);
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::consistent_second_invariant_of_deviatoric_tensor(shear_strain_rate), 0.));
 
               const double adiabatic_temperature = this->get_adiabatic_conditions().is_initialized()
                                                    ?

source/material_model/reaction_model/grain_size_evolution.cc

@@ -226,8 +226,8 @@ namespace aspect
                                           std::pow(roughness_to_grain_size, m);
 
               // 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 SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::consistent_deviator(in.strain_rate[i]);
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::consistent_second_invariant_of_deviatoric_tensor(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::consistent_second_invariant_of_deviatoric_tensor(Utilities::Tensors::consistent_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::consistent_second_invariant_of_deviatoric_tensor(Utilities::Tensors::consistent_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::consistent_second_invariant_of_deviatoric_tensor(Utilities::Tensors::consistent_deviator(strain_rate)), 0.)),
                                         min_strain_rate);
         const double log_edot_ii = std::log(edot_ii);