change the names of the modified tensor functions

YiminJin · YiminJin · commit 447f09f09628 · 2025-06-17T15:09:45.000-07:00
diff --git a/doc/modules/changes/20250617_yiminjin b/doc/modules/changes/20250617_yiminjin
@@ -1,5 +1,5 @@
-Changed: The dealii functions deviator() and second_invariant() are replaced by
-Utilities::Tensors::deviator() and Utilities::Tensors::deviatoric_tensor_inv2(),
-which are consistent with the plane strain assumption in 2D.
+Changed: The deviator and second invariant of symmetric tensors are modified to
+be consistent with the plane strain assumption in 2D. It changes the outputs of
+compressible material models that are dependent on the deviatoric strain rate.
 <br>
 (Yimin Jin, 2025/06/17)
diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h
@@ -1409,27 +1409,27 @@ namespace aspect
       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.
+       * 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.
        */
       template <int dim>
       SymmetricTensor<2,dim>
-      deviator(const SymmetricTensor<2,dim> &input);
+      plane_strain_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
+       * 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.
        */
       template <int dim>
       double
-      deviatoric_tensor_inv2(const SymmetricTensor<2,dim> &input);
+      plane_strain_second_invariant(const SymmetricTensor<2,dim> &input);
     }
 
   }
diff --git a/source/heating_model/shear_heating.cc b/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(Utilities::Tensors::deviatoric_tensor_inv2(deviatoric_strain_rate)));
+              const double deviatoric_stress = 2 * material_model_outputs.viscosities[q] * std::sqrt(std::fabs(Utilities::Tensors::plane_strain_second_invariant(deviatoric_strain_rate)));
               double scaling_factor = 1.0;
               if (deviatoric_stress > 0.0)
                 scaling_factor = std::min(yield_stress / deviatoric_stress, 1.0);
diff --git a/source/material_model/drucker_prager.cc b/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 = Utilities::Tensors::deviator(in.strain_rate[i]);
+              const SymmetricTensor<2,dim> strain_rate_deviator = Utilities::Tensors::plane_strain_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
@@ -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(Utilities::Tensors::deviatoric_tensor_inv2(strain_rate_deviator))
+                                             std::fabs(Utilities::Tensors::plane_strain_second_invariant(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(Utilities::Tensors::deviatoric_tensor_inv2(strain_rate_deviator))));
+                                              std::fabs(Utilities::Tensors::plane_strain_second_invariant(strain_rate_deviator))));
 
               const double strain_rate_effective = edot_ii_strict;
 
diff --git a/source/material_model/entropy_model.cc b/source/material_model/entropy_model.cc
@@ -253,7 +253,7 @@ namespace aspect
                                                                                                       drucker_prager_parameters);
                     }
 
-                  const double strain_rate_effective = std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])));
+                  const double strain_rate_effective = std::fabs(Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(in.strain_rate[i])));
 
                   if (std::sqrt(strain_rate_effective) >= std::numeric_limits<double>::min())
                     {
diff --git a/source/material_model/grain_size.cc b/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 = Utilities::Tensors::deviator(strain_rate);
-      const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(shear_strain_rate), 0.));
+      const SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::plane_strain_deviator(strain_rate);
+      const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(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(-Utilities::Tensors::deviatoric_tensor_inv2(dislocation_strain_rate), 0.));
+          const double dislocation_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(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 = 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 SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::plane_strain_deviator(in.strain_rate[i]);
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(shear_strain_rate), 0.));
 
               const double adiabatic_temperature = this->get_adiabatic_conditions().is_initialized()
                                                    ?
diff --git a/source/material_model/reaction_model/grain_size_evolution.cc b/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 = 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 SymmetricTensor<2,dim> shear_strain_rate = Utilities::Tensors::plane_strain_deviator(in.strain_rate[i]);
+              const double second_strain_rate_invariant = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(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);
diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/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(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_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(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate)), 0.)),
                                         minimum_strain_rate);
         const double log_edot_ii = std::log(edot_ii);
 
diff --git a/source/material_model/rheology/diffusion_dislocation.cc b/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(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.)),
+        const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate)), 0.)),
                                         min_strain_rate);
         const double log_edot_ii = std::log(edot_ii);
 
diff --git a/source/material_model/rheology/elasticity.cc b/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 = Utilities::Tensors::deviator(in.strain_rate[i]);
+                const SymmetricTensor<2, dim> deviatoric_strain_rate = Utilities::Tensors::plane_strain_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 = Utilities::Tensors::deviator(in.strain_rate[i]);
+            const SymmetricTensor<2, dim> deviatoric_strain_rate = Utilities::Tensors::plane_strain_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 * Utilities::Tensors::deviator(in.strain_rate[i])
+                stress_t = 2. * effective_creep_viscosity * Utilities::Tensors::plane_strain_deviator(in.strain_rate[i])
                            + effective_creep_viscosity / elastic_viscosity * stress_0_t
                            + (1. - timestep_ratio) * (1. - effective_creep_viscosity / elastic_viscosity) * stress_old;
 
diff --git a/source/material_model/rheology/strain_dependent.cc b/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(Utilities::Tensors::deviatoric_tensor_inv2(L));
+              const double strain_ii = std::fabs(Utilities::Tensors::plane_strain_second_invariant(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(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])), 0.)),
+            const double edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(in.strain_rate[i])), 0.)),
                                             min_strain_rate);
             double delta_e_ii = edot_ii*this->get_timestep();
 
diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc
@@ -153,7 +153,7 @@ namespace aspect
           edot_ii = ref_strain_rate;
         else
           // Calculate the square root of the second moment invariant for the deviatoric strain rate tensor.
-          edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])), 0.)),
+          edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(in.strain_rate[i])), 0.)),
                              min_strain_rate);
 
         // Calculate viscosities for each of the individual compositional phases
@@ -307,13 +307,13 @@ namespace aspect
                                       "not filled by the caller."));
 
                     const SymmetricTensor<2, dim> effective_strain_rate =
-                      elastic_rheology.calculate_viscoelastic_strain_rate(Utilities::Tensors::deviator(in.strain_rate[i]),
+                      elastic_rheology.calculate_viscoelastic_strain_rate(Utilities::Tensors::plane_strain_deviator(in.strain_rate[i]),
                                                                           stress_0_advected,
                                                                           stress_old,
                                                                           non_yielding_viscosity,
                                                                           elastic_shear_modulus);
 
-                    effective_edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(effective_strain_rate), 0.)),
+                    effective_edot_ii = std::max(std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(effective_strain_rate), 0.)),
                                                  min_strain_rate);
                   }
               }
diff --git a/source/particle/property/integrated_strain_invariant.cc b/source/particle/property/integrated_strain_invariant.cc
@@ -62,7 +62,7 @@ namespace aspect
             const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u);
 
             // Calculate strain rate second invariant
-            const double edot_ii = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.));
+            const double edot_ii = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate)), 0.));
 
             // New strain is the old strain plus dt*edot_ii
             const double new_strain = old_strain + dt*edot_ii;
diff --git a/source/particle/property/viscoplastic_strain_invariants.cc b/source/particle/property/viscoplastic_strain_invariants.cc
@@ -137,7 +137,7 @@ namespace aspect
             const auto data = particle.get_properties();
 
             // Calculate strain rate second invariant
-            const double edot_ii = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(material_inputs.strain_rate[0])), 0.));
+            const double edot_ii = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(material_inputs.strain_rate[0])), 0.));
 
             // Calculate strain invariant magnitude over the last time step
             const double strain_update = dt*edot_ii;
diff --git a/source/postprocess/boundary_strain_rate_residual_statistics.cc b/source/postprocess/boundary_strain_rate_residual_statistics.cc
@@ -126,7 +126,7 @@ namespace aspect
                       const Point<dim> point_at_surface = fe_face_values.quadrature_point(q);
 
                       // Extract the strain rate invariant output here.
-                      const double e_ii = std::sqrt(std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate[q]))));
+                      const double e_ii = std::sqrt(std::fabs(Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate[q]))));
 
                       // Extract data strain rate invartiant.
                       double e_ii_data = get_data_surface_strain_rate(point_at_surface);
diff --git a/source/postprocess/visualization/boundary_strain_rate_residual.cc b/source/postprocess/visualization/boundary_strain_rate_residual.cc
@@ -100,7 +100,7 @@ namespace aspect
               // Only compute residual for doubles. This condition checks for nan values
               if (data_surface_strain_rate < 1e300 || !std::isnan(data_surface_strain_rate))
                 computed_quantities[q](0) = data_surface_strain_rate -
-                                            std::sqrt(std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)))) * unit_scaling_factor;
+                                            std::sqrt(std::fabs(Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate)))) * unit_scaling_factor;
               else
                 continue;
 
diff --git a/source/postprocess/visualization/strain_rate.cc b/source/postprocess/visualization/strain_rate.cc
@@ -61,7 +61,7 @@ namespace aspect
               grad_u[d] = input_data.solution_gradients[q][d];
 
             const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u);
-            computed_quantities[q](0) = std::sqrt(std::max(-Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(strain_rate)), 0.));
+            computed_quantities[q](0) = std::sqrt(std::max(-Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(strain_rate)), 0.));
           }
 
         // average the values if requested
diff --git a/source/postprocess/visualization/stress_residual.cc b/source/postprocess/visualization/stress_residual.cc
diff --git a/source/postprocess/visualization/stress_second_invariant.cc b/source/postprocess/visualization/stress_second_invariant.cc
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc
diff --git a/source/utilities.cc b/source/utilities.cc
diff --git a/tests/melt_visco_plastic.cc b/tests/melt_visco_plastic.cc
diff --git a/tests/spiegelman_fail_test.cc b/tests/spiegelman_fail_test.cc

Original file line number	Diff line number	Diff line change
`@@ -253,7 +253,7 @@ namespace aspect`
`253`	`253`	`drucker_prager_parameters);`
`254`	`254`	`}`
`255`	`255`
`256`		`- const double strain_rate_effective = std::fabs(Utilities::Tensors::deviatoric_tensor_inv2(Utilities::Tensors::deviator(in.strain_rate[i])));`
	`256`	`+ const double strain_rate_effective = std::fabs(Utilities::Tensors::plane_strain_second_invariant(Utilities::Tensors::plane_strain_deviator(in.strain_rate[i])));`
`257`	`257`
`258`	`258`	`if (std::sqrt(strain_rate_effective) >= std::numeric_limits<double>::min())`
`259`	`259`	`{`