Fix elastic stress used in principal_stress postprocessor #5480

gassmoeller · 2023-11-07T13:42:46Z

Similar to #4375, but for the principal_stress postprocessor.

@rfildes: Could you check if this change also fixes the bug you observed?

@bobmyhill or @anne-glerum: Could one of you review this PR?

bangerth

The postprocessor is documented as follows:

                                                  "A visualization output object that outputs the "
                                                  "principal stress values and directions, i.e., the "
                                                  "eigenvalues and eigenvectors of the stress tensor. "
                                                  "The postprocessor can either operate on the full "
                                                  "stress tensor or only on the deviatoric stress tensor, "
                                                  "depending on what run-time parameters are set."

I think it would be useful to define what "principal stress" actually means here, so that we can reason what the postprocessor should actually do. As it stands, it's hard to tell whether what you do is correct or not because we haven't defined what this postprocessor should actually do.

bangerth · 2023-11-08T16:14:52Z

source/postprocess/visualization/principal_stress.cc

            if (this->get_parameters().enable_elasticity == true)
              {
-                stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
-                stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
-                stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
+                stress[0][0] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
+                stress[1][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
+                stress[0][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

                if (dim == 3)
                  {
-                    stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
-                    stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
-                    stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
+                    stress[2][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
+                    stress[0][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
+                    stress[1][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
                  }
              }
+            else
+              {
+                const SymmetricTensor<2,dim> deviatoric_strain_rate
+                  = (this->get_material_model().is_compressible()
+                     ?
+                     in.strain_rate[q] - 1./3. * trace(in.strain_rate[q]) * unit_symmetric_tensor<dim>()
+                     :
+                     in.strain_rate[q]);
+
+                // Compressive stress is positive in geoscience applications
+                stress = -2. * out.viscosities[q] * deviatoric_strain_rate;
+              }


Would you mind switching the if and else branches here? Most readers will approach this from the non-elasticity case in mind, and it would be nice to see the non-elastic case first.

I think also having the commented note about sign convention (the need to convert from the convention used inside of the code to that used in geology and in post processing here) moved above both the viscous and viscoelastic cases would make it clearer what is being done

rfildes · 2023-11-09T18:18:20Z

source/postprocess/visualization/principal_stress.cc

+                stress[0][0] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
+                stress[1][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
+                stress[0][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];


I think the sign should be changed here similarly in how it is changed for the viscous case in line 144 to convert from the convention used inside of the code to the geologic convention of compression as positive

Suggested change

stress[0][0] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];

stress[1][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];

stress[0][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

stress[0][0] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];

stress[1][1] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];

stress[0][1] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

rfildes · 2023-11-09T18:18:27Z

source/postprocess/visualization/principal_stress.cc

+                    stress[2][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
+                    stress[0][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
+                    stress[1][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];


Suggested change

stress[2][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];

stress[0][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];

stress[1][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];

stress[2][2] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];

stress[0][2] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];

stress[1][2] = - in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];

rfildes · 2023-11-09T18:33:22Z

@gassmoeller I have tested this and found it does not fix the stress orientation issue in the hanging beam test ( tests and description on this thread: https://community.geodynamics.org/t/principal-stress-orientations-in-viscous-vs-visco-elastic-models/3161/5). I made some comments to your code changes above where I think the remaining issue is.

gassmoeller · 2023-11-09T20:09:30Z

All good points, thanks for the review. I have corrected the signs and added explanation to the documentation of the postprocessor. The problem with the sign is really that it is a convention, and we just need to be consistent to use the same convention.
@rfildes do you want to contribute your test cases yourself in a separate PR (this would give you more credit for it), or do you want me to incorporate them into this PR?

rfildes · 2023-11-10T01:16:28Z

@gassmoeller You can incorporate those test cases with this PR. I think that makes sense so it is all together! (I also was having some technical problems on my end and don't want to hold up the merging)

gassmoeller · 2023-11-16T18:43:59Z

Ok, I added your tests to this PR and made sure the results are similar to the ones you posted in the forum.

This fix should be ready for a final review and merge.

bangerth · 2023-11-19T12:46:05Z

One of the testers is stuck, but since the others are all succeeding, I'll take the liberty to merge.

gassmoeller force-pushed the fix_principal_stress_elastic branch 2 times, most recently from ce7921a to 54f3af7 Compare November 8, 2023 15:32

bangerth reviewed Nov 8, 2023

View reviewed changes

rfildes reviewed Nov 9, 2023

View reviewed changes

gassmoeller force-pushed the fix_principal_stress_elastic branch from 54f3af7 to fe7d12d Compare November 9, 2023 20:02

gassmoeller force-pushed the fix_principal_stress_elastic branch from ec1f3bb to 9efbd64 Compare November 16, 2023 16:49

Fix stress used in principal_stress.cc

b8442ef

gassmoeller force-pushed the fix_principal_stress_elastic branch from 9efbd64 to b8442ef Compare November 16, 2023 16:51

bangerth approved these changes Nov 19, 2023

View reviewed changes

bangerth merged commit 299a645 into geodynamics:main Nov 19, 2023
7 of 8 checks passed

gassmoeller deleted the fix_principal_stress_elastic branch November 27, 2023 16:21

bobmyhill mentioned this pull request Nov 28, 2023

fixed stress visualization for elastic materials #4375

Merged

4 tasks

gassmoeller restored the fix_principal_stress_elastic branch December 8, 2023 18:18

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix elastic stress used in principal_stress postprocessor #5480

Fix elastic stress used in principal_stress postprocessor #5480

gassmoeller commented Nov 7, 2023

bangerth left a comment

bangerth Nov 8, 2023

rfildes Nov 9, 2023

rfildes Nov 9, 2023 •

edited

rfildes Nov 9, 2023

rfildes commented Nov 9, 2023

gassmoeller commented Nov 9, 2023

rfildes commented Nov 10, 2023

gassmoeller commented Nov 16, 2023

bangerth commented Nov 19, 2023

Fix elastic stress used in principal_stress postprocessor #5480

Fix elastic stress used in principal_stress postprocessor #5480

Conversation

gassmoeller commented Nov 7, 2023

bangerth left a comment

Choose a reason for hiding this comment

bangerth Nov 8, 2023

Choose a reason for hiding this comment

rfildes Nov 9, 2023

Choose a reason for hiding this comment

rfildes Nov 9, 2023 • edited

Choose a reason for hiding this comment

rfildes Nov 9, 2023

Choose a reason for hiding this comment

rfildes commented Nov 9, 2023

gassmoeller commented Nov 9, 2023

rfildes commented Nov 10, 2023

gassmoeller commented Nov 16, 2023

bangerth commented Nov 19, 2023

rfildes Nov 9, 2023 •

edited