added shear stress postprocessor

geodynamics · Nov 10, 2023 · 3391ba8 · 3391ba8
1 parent 3ea11c6
commit 3391ba8
Show file tree

Hide file tree

Showing 6 changed files with 45,041 additions and 45,031 deletions.
diff --git a/doc/modules/changes/20231110_bobmyhill b/doc/modules/changes/20231110_bobmyhill
@@ -1,4 +1,4 @@
-Fixed: The stress visualization postprocessor now outputs
-the correct stresses when elasticity is enabled.
+Fixed: The stress and shear stress visualization postprocessors
+now output the correct stresses when elasticity is enabled.
 <br>
-(Bob Myhill, 2023/11/10)
+(Bob Myhill, Rebecca Fildes, 2023/11/10)
diff --git a/source/postprocess/visualization/shear_stress.cc b/source/postprocess/visualization/shear_stress.cc
@@ -67,33 +67,43 @@ namespace aspect
         // ...and use it to compute the stresses
         for (unsigned int q=0; q<n_quadrature_points; ++q)
           {
-            const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
-            const SymmetricTensor<2,dim> deviatoric_strain_rate
-              = (this->get_material_model().is_compressible()
-                 ?
-                 strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
-                 :
-                 strain_rate);
-
-            const double eta = out.viscosities[q];
-
-            // Compressive stress is positive in geoscience applications
-            SymmetricTensor<2,dim> shear_stress = -2.*eta*deviatoric_strain_rate;
-
-            // Add elastic stresses if existent
+            // Compressive stress is negative by the sign convention
+            // used by the engineering community, and internally by ASPECT.
+            // Here, we change the sign of the stress to match the
+            // sign convention used by the geoscience community,
+            // which is also the sign convention for stress expected
+            // in ASPECT parameter files.
+            SymmetricTensor<2,dim> shear_stress;
+
+            // If elasticity is enabled, the deviatoric stress is stored
+            // in compositional fields, otherwise the deviatoric stress
+            // can be obtained from the viscosity and strain rate.
             if (this->get_parameters().enable_elasticity == true)
               {
-                shear_stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
-                shear_stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
-                shear_stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
+                shear_stress[0][0] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
+                shear_stress[1][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
+                shear_stress[0][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
 
                 if (dim == 3)
                   {
-                    shear_stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
-                    shear_stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
-                    shear_stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
+                    shear_stress[2][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
+                    shear_stress[0][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
+                    shear_stress[1][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
                   }
               }
+            else
+              {
+                const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
+                const SymmetricTensor<2,dim> deviatoric_strain_rate
+                  = (this->get_material_model().is_compressible()
+                     ?
+                     strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
+                     :
+                     strain_rate);
+
+                const double eta = out.viscosities[q];
+                shear_stress -= 2.*eta*deviatoric_strain_rate;
+              }
 
             for (unsigned int d=0; d<dim; ++d)
               for (unsigned int e=0; e<dim; ++e)

diff --git a/source/postprocess/visualization/stress.cc b/source/postprocess/visualization/stress.cc
@@ -102,7 +102,7 @@ namespace aspect
                      strain_rate);
 
                 const double eta = out.viscosities[q];
-                stress += -2.*eta*deviatoric_strain_rate;
+                stress -= 2.*eta*deviatoric_strain_rate;
               }
 
             for (unsigned int d=0; d<dim; ++d)

diff --git a/tests/viscoelastic_stress_build-up.prm b/tests/viscoelastic_stress_build-up.prm
@@ -77,7 +77,7 @@ end
 subsection Postprocess
   subsection Visualization
     set Time between graphical output = 1e3
-    set List of output variables = stress second invariant
+    set List of output variables = shear stress, stress second invariant
     set Output format = gnuplot
   end
 end