Add cpo elastic tensor decomposition plugin #5116

MFraters · 2023-03-21T01:09:15Z

Last major part of #3885. Adds a property plugin which can use the elastic tensor computed in #4987 to compute actual useful values from it. It requires that plugin to be present, so it is currently rebased on that.

Since this is the last major component of the CPO plugin set, I am happy to include the integration tests within this pull request, or do that in a follow up pull request. I will include the changelog entry in the pull request where the integration tests are added, so if that is here, that will be here as well.

This pull request is ready for review, but can't be merged until #4987 is merged.

bobmyhill · 2023-07-08T06:37:46Z

Hi @MFraters, mostly typos here, but I do wonder whether the vast majority of this shouldn't be in an TensorConversion namespace in utilities rather than a separate plugin. Most of the functions are quite general, and they would probably be useful for anisotropic viscosity too, which also describes a linear relationship between two symmetric second order tensors.

You should extend the docstring(s) introducing the conversions, the symmetries behind them, and citing the equations from a relevant study for each (probably Mainprice). Expecially for the conversions to and from the 1x21 vector, which I don't believe is a standard form.

P.S.: Most of my comments here mirror @gassmoeller's from #4987. Good to double check his comments as well as mine to make sure you've covered everything.

include/aspect/particle/property/cpo_elastic_tensor.h

bobmyhill · 2023-07-07T18:05:46Z

include/aspect/particle/property/cpo_elastic_tensor.h

+    {
+
+      /**
+       * Computes the elastic tensor $C_{ijkl} based on the cpo in both olivine


This is documentation for a class, not a function. Can you reword the docstring to clarify the purpose of the class?

bobmyhill · 2023-07-07T18:19:01Z

include/aspect/particle/property/cpo_elastic_tensor.h

+           * grain size, $N_s$ is the number of grains and $C_{ijkl}$ is the elastic
+           * tensor. This elastic tensor is computed by the equation above in
+           * Mainprice (1990): $C_{ijkl} = R_{ip} R_{jg} R_{kr} R_{is} C_{pgrs}$, where
+           * R_{ij} is the cpo orientation matrix.


Can you mention that you return the Voigt form?

bobmyhill · 2023-07-07T19:15:12Z

include/aspect/particle/property/cpo_elastic_tensor.h

+                             SymmetricTensor<2,6> &elastic_tensor);
+
+          /**
+           * Rotate a 3D 4th order tensor with an other 3D 4th


Suggested change

* Rotate a 3D 4th order tensor with an other 3D 4th

* Rotate a 3x3x3x3 elastic tensor with a 3x3 tensor

bobmyhill · 2023-07-07T19:15:57Z

include/aspect/particle/property/cpo_elastic_tensor.h

+
+
+          /**
+           * Rotate a 6x6 voigt matrix with an other 3D 4th


Suggested change

* Rotate a 6x6 voigt matrix with an other 3D 4th

* Rotate a 4th order elastic tensor in Voigt form (6x6) with a 3x3 rotation tensor.

include/aspect/particle/property/elastic_tensor_decomposition.h

bobmyhill · 2023-07-08T05:37:42Z

source/particle/property/cpo_elastic_tensor.cc

+                                                           const ArrayView<double> &data) const
+      {
+        SymmetricTensor<2,6,double> Sav;
+        const SymmetricTensor<2,6> *stiffness_matrix = &stiffness_matrix_olivine;


there's an odd mix of hard coding here and generality in the next line

bobmyhill · 2023-07-08T06:15:48Z

source/particle/property/cpo_elastic_tensor.cc

+      template <int dim>
+      CpoElasticTensor<dim>::CpoElasticTensor ()
+      {
+        permutation_operator_3d[0][1][2]  = 1;


This is the Levi-Civita tensor and is already used in cpo_bingham_average. Better to move somewhere more general.

bobmyhill · 2023-07-08T06:27:17Z

source/particle/property/cpo_elastic_tensor.cc

+          input[0][0],           // 0  // 1
+          input[1][1],           // 1  // 2
+          input[2][2],           // 2  // 3
+          sqrt(2)*input[1][2],   // 3  // 4


Please use the std namespace for sqrt and similar math functions.

bobmyhill · 2023-07-08T06:28:59Z

unit_tests/crystal_preferred_orientation.cc

+      double alpha = radians;
+      double beta = radians;
+      double gamma = radians;
+      rotation_tensor[0][0] = cos(alpha) * cos(beta);


Please use the std namespace for sin, cos and similar math functions.

MFraters

@bobmyhill Thanks for your review! I have addressed your comments related to the files for this pull request.

MFraters · 2023-08-03T21:17:25Z

include/aspect/particle/property/elastic_tensor_decomposition.h

+
+          /**
+           * Return which data has to be provided to update the property.
+           * The integrated strains needs the gradients of the velocity.


I removed the function, since it just uses the default

MFraters · 2023-08-03T21:26:27Z

include/aspect/particle/property/elastic_tensor_decomposition.h

+          indexed_even_permutation(const unsigned short int index);
+
+          /**
+           * Computes the voigt stiffness tensor from the elastic tensor


Why do you expect it to be 6x6?

include/aspect/particle/property/elastic_tensor_decomposition.h

gassmoeller

First batch of comments. I can take another look after you moved the matrices.

include/aspect/particle/property/elastic_tensor_decomposition.h

gassmoeller · 2023-10-08T00:20:51Z

source/particle/property/elastic_tensor_decomposition.cc

+                if (std::abs(dv_dot_product) >= 1.0)
+                  dv_dot_product = std::copysign(1.0,dv_dot_product);
+                // compute the angle bewteen the vectors and account for that vectors in the oposit
+                // direction are the same. So limit them between 0 and 90 degrees such that it


this first sentence doesnt really make sense. please rework.

I hope this is more clear.

source/particle/property/elastic_tensor_decomposition.cc

MFraters

Thanks for your review and sorry that it took so long to get around to this. I have addressed all your comments.

MFraters · 2024-03-16T03:08:21Z

source/particle/property/elastic_tensor_decomposition.cc

+        projection_matrix_tric_to_mono[7][7] = 1.0;
+        projection_matrix_tric_to_mono[8][8] = 1.0;
+        projection_matrix_tric_to_mono[11][11] = 1.0;
+        projection_matrix_tric_to_mono[14][14] = 1.0;
+        projection_matrix_tric_to_mono[17][17] = 1.0;
+        projection_matrix_tric_to_mono[20][20] = 1.0;
+
+
+        // projection_matrix_mono_to_ortho;
+        projection_matrix_mono_to_ortho[0][0] = 1.0;
+        projection_matrix_mono_to_ortho[1][1] = 1.0;
+        projection_matrix_mono_to_ortho[2][2] = 1.0;
+        projection_matrix_mono_to_ortho[3][3] = 1.0;
+        projection_matrix_mono_to_ortho[4][4] = 1.0;
+        projection_matrix_mono_to_ortho[5][5] = 1.0;
+        projection_matrix_mono_to_ortho[6][6] = 1.0;
+        projection_matrix_mono_to_ortho[7][7] = 1.0;
+        projection_matrix_mono_to_ortho[8][8] = 1.0;
+
+        // projection_matrix_ortho_to_tetra;
+        projection_matrix_ortho_to_tetra[0][0] = 0.5;
+        projection_matrix_ortho_to_tetra[0][1] = 0.5;
+        projection_matrix_ortho_to_tetra[1][1] = 0.5;
+        projection_matrix_ortho_to_tetra[2][2] = 1.0;
+        projection_matrix_ortho_to_tetra[3][3] = 0.5;
+        projection_matrix_ortho_to_tetra[3][4] = 0.5;
+        projection_matrix_ortho_to_tetra[4][4] = 0.5;
+        projection_matrix_ortho_to_tetra[5][5] = 1.0;
+        projection_matrix_ortho_to_tetra[6][6] = 0.5;
+        projection_matrix_ortho_to_tetra[6][7] = 0.5;
+        projection_matrix_ortho_to_tetra[7][7] = 0.5;
+        projection_matrix_ortho_to_tetra[8][8] = 1.0;
+
+        // projection_matrix_tetra_to_hexa;
+        projection_matrix_tetra_to_hexa[0][0] = 3./8.;
+        projection_matrix_tetra_to_hexa[0][1] = 3./8.;
+        projection_matrix_tetra_to_hexa[1][1] = 3./8.;
+        projection_matrix_tetra_to_hexa[2][2] = 1.0;
+        projection_matrix_tetra_to_hexa[3][3] = 0.5;
+        projection_matrix_tetra_to_hexa[3][4] = 0.5;
+        projection_matrix_tetra_to_hexa[4][4] = 0.5;
+        projection_matrix_tetra_to_hexa[5][5] = 3./4.;
+        projection_matrix_tetra_to_hexa[6][6] = 0.5;
+        projection_matrix_tetra_to_hexa[6][7] = 0.5;
+        projection_matrix_tetra_to_hexa[7][7] = 0.5;
+        projection_matrix_tetra_to_hexa[8][8] = 0.5;
+        projection_matrix_tetra_to_hexa[5][0] = 1./(4.*std::sqrt(2.0));
+        projection_matrix_tetra_to_hexa[5][1] = 1./(4.*std::sqrt(2.0));
+        projection_matrix_tetra_to_hexa[8][0] = 0.25;
+        projection_matrix_tetra_to_hexa[8][1] = 0.25;
+        projection_matrix_tetra_to_hexa[8][5] = -1/(2*std::sqrt(2.0));
+
+
+        // projection_matrix_hexa_to_iso;
+        projection_matrix_hexa_to_iso[0][0] = 3./15.;
+        projection_matrix_hexa_to_iso[0][1] = 3./15.;
+        projection_matrix_hexa_to_iso[0][2] = 3./15.;
+        projection_matrix_hexa_to_iso[1][1] = 3./15.;
+        projection_matrix_hexa_to_iso[1][2] = 3./15.;
+        projection_matrix_hexa_to_iso[2][2] = 3./15.;
+        projection_matrix_hexa_to_iso[3][3] = 4./15.;
+        projection_matrix_hexa_to_iso[3][4] = 4./15.;
+        projection_matrix_hexa_to_iso[3][5] = 4./15.;
+        projection_matrix_hexa_to_iso[4][4] = 4./15.;
+        projection_matrix_hexa_to_iso[4][5] = 4./15.;
+        projection_matrix_hexa_to_iso[5][5] = 4./15.;
+        projection_matrix_hexa_to_iso[6][6] = 1./5.;
+        projection_matrix_hexa_to_iso[6][7] = 1./5.;
+        projection_matrix_hexa_to_iso[6][8] = 1./5.;
+        projection_matrix_hexa_to_iso[7][7] = 1./5.;
+        projection_matrix_hexa_to_iso[7][8] = 1./5.;
+        projection_matrix_hexa_to_iso[8][8] = 1./5.;
+
+        projection_matrix_hexa_to_iso[0][3] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[0][4] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[0][5] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[1][3] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[1][4] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[1][5] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[2][3] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[2][4] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[2][5] = std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[0][6] = 2./15.;
+        projection_matrix_hexa_to_iso[0][7] = 2./15.;
+        projection_matrix_hexa_to_iso[0][8] = 2./15.;
+        projection_matrix_hexa_to_iso[1][6] = 2./15.;
+        projection_matrix_hexa_to_iso[1][7] = 2./15.;
+        projection_matrix_hexa_to_iso[1][8] = 2./15.;
+        projection_matrix_hexa_to_iso[2][6] = 2./15.;
+        projection_matrix_hexa_to_iso[2][7] = 2./15.;
+        projection_matrix_hexa_to_iso[2][8] = 2./15.;
+        projection_matrix_hexa_to_iso[3][6] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[3][7] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[3][8] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[4][6] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[4][7] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[4][8] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[5][6] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[5][7] = -std::sqrt(2.0)/15.;
+        projection_matrix_hexa_to_iso[5][8] = -std::sqrt(2.0)/15.;


I can't make them constexpr, but I can make them const, so that is what I did. Feel free to see if you can make it work for constexpr, but I wasn't able to

MFraters · 2024-03-16T04:02:29Z

source/particle/property/elastic_tensor_decomposition.cc

+                if (std::abs(dv_dot_product) >= 1.0)
+                  dv_dot_product = std::copysign(1.0,dv_dot_product);
+                // compute the angle bewteen the vectors and account for that vectors in the oposit
+                // direction are the same. So limit them between 0 and 90 degrees such that it


I hope this is more clear.

source/particle/property/elastic_tensor_decomposition.cc

include/aspect/particle/property/elastic_tensor_decomposition.h

MFraters · 2024-05-29T16:30:40Z

@gassmoeller Could you have another look at this?

gassmoeller

The structure looks good now, I only have a bunch of small comments and requests for better documentation. Can you also add a changelog entry?

Also: We dont have many tests in the tests/ directory for the rest of CPO implementation. I vaguely remember that you said you need this PR to be able to write tests. Can you start adding tests to this PR? Otherwise we will forget about it in the future.

gassmoeller · 2024-05-30T03:20:16Z

include/aspect/particle/property/elastic_tensor_decomposition.h

+
+      namespace Utility
+      {
+


Have a look and be consistent about the number of empty lines when opening namespaces. I would suggest no empty lines in line 39, and only one in line 32-33.

include/aspect/particle/property/elastic_tensor_decomposition.h

gassmoeller · 2024-05-30T20:00:16Z

source/particle/property/elastic_tensor_decomposition.cc

+            // The following line would do the same as the lines below, but it is slow. It has therefore been
+            // replaced by the lines below.
+            //auto ortho_and_higher_vector = projection_matrix_mono_to_ortho*mono_and_higher_vector;
+            dealii::Tensor<1,9>  mono_and_higher_vector_cropped;


What does cropped mean?

gassmoeller · 2024-05-30T20:05:26Z

source/particle/property/elastic_tensor_decomposition.cc

+            ortho_and_higher_vector[6] = mono_and_higher_vector[6];
+            ortho_and_higher_vector[7] = mono_and_higher_vector[7];
+            ortho_and_higher_vector[8] = mono_and_higher_vector[8];
+            auto mono_vector = mono_and_higher_vector_cropped-ortho_and_higher_vector;


can you change the variable names to always contain the full name of the crystal system? mono is a very ambiguous name, monoclinal is much more specific.

gassmoeller · 2024-05-30T20:10:35Z

source/particle/property/elastic_tensor_decomposition.cc

+
+      template <int dim>
+      void
+      ElasticTensorDecomposition<dim>::declare_parameters (ParameterHandler &)


if declare_parameters and parse_parameters are empty remove them

source/particle/property/elastic_tensor_decomposition.cc

gassmoeller · 2024-05-30T20:11:16Z

source/particle/property/elastic_tensor_decomposition.cc

+                                        "elastic tensor decomposition",
+                                        "A plugin in which decomposes the elastic tensor "
+                                        "into different approximations and provides the "
+                                        "eigenvectors of the tensor.")


Is it worth describing in a bit more details what the different approximations are?

MFraters · 2024-06-04T20:57:36Z

@gassmoeller I addressed all your comments, except for adding a test. I will do that next, but the rest is ready for review again.

gassmoeller · 2024-06-04T22:16:14Z

There are a number of my comments that are marked as resolved, but have not been fixed in the code (it is a bad habit of Github to mark comments as resolved if you rebase and push, I think this is why they dont show up). Can you check https://github.com/geodynamics/aspect/pull/5116/files and check all resolved comments if they have been addressed?

Co-authored-by: Bob Myhill <[email protected]>

Co-authored-by: Rene Gassmoeller <[email protected]>

MFraters · 2024-06-04T22:29:10Z

Sorry, I think I forgot to pull after applying your comments and before I added the changes locally and then force pushed. I fixed it now manually.

gassmoeller

Alright, I think this is ready to be included. Last major part of the CPO implementation done. Thanks again @MFraters for undertaking this massive project! Looking forward to seeing more applications using CPO. 🎉

As a follow-up PR: Think about removing the warnings that CPO is experimental. Add more integration tests, and a changelog entry?

gassmoeller · 2024-06-05T03:23:43Z

/rebuild

MFraters · 2024-06-05T04:43:57Z

Thanks for all the reviewing, it is really appreciated and made the code a whole lot better! 🎉

bobmyhill reviewed Jul 8, 2023

View reviewed changes

MFraters mentioned this pull request Jul 13, 2023

Add elastic utilities #5284

Merged

MFraters mentioned this pull request Jul 13, 2023

Tracking issue for adding CPO/LPO tracking and development to aspect #3885

Open

MFraters force-pushed the add_cpo_elastic_tensor_decomposition_plugin branch 3 times, most recently from 0dd656d to f88113a Compare August 4, 2023 19:12

MFraters commented Aug 4, 2023

View reviewed changes

MFraters force-pushed the add_cpo_elastic_tensor_decomposition_plugin branch from f88113a to f8a4bcd Compare August 12, 2023 19:31

MFraters mentioned this pull request Aug 14, 2023

CPO code needs more tests #5361

Closed

gassmoeller reviewed Oct 8, 2023

View reviewed changes

MFraters commented Mar 16, 2024

View reviewed changes

MFraters force-pushed the add_cpo_elastic_tensor_decomposition_plugin branch from e946b52 to 9fd6dc1 Compare March 16, 2024 04:13

MFraters changed the title ~~[ready for review, do no merge] Add cpo elastic tensor decomposition plugin~~ Add cpo elastic tensor decomposition plugin Mar 16, 2024

gassmoeller self-assigned this Apr 11, 2024

gassmoeller reviewed May 30, 2024

View reviewed changes

Wang-yijun mentioned this pull request May 31, 2024

Change plugin name in cpo elastic tensor and use olivine stiffness matrix in initialization #5690

Merged

MFraters force-pushed the add_cpo_elastic_tensor_decomposition_plugin branch 2 times, most recently from b0e31f6 to 55be865 Compare June 4, 2024 20:57

MFraters and others added 6 commits June 4, 2024 16:27

Add particle property elastic tensor decomposition plugin.

afdfc97

add unit test for elastic tensor decomposition particle property plugin.

f3fdce4

Apply suggestions from code review from Bob

e390510

Co-authored-by: Bob Myhill <[email protected]>

Further adress Bob's comments.

71dac2b

Apply suggestions from code review Rene

d9c4b99

Co-authored-by: Rene Gassmoeller <[email protected]>

Adress comments Rene

7368f89

MFraters force-pushed the add_cpo_elastic_tensor_decomposition_plugin branch from 55be865 to 7368f89 Compare June 4, 2024 22:27

gassmoeller approved these changes Jun 5, 2024

View reviewed changes

gassmoeller added ready to merge ready to test labels Jun 5, 2024

gassmoeller merged commit 3c22a3a into geodynamics:main Jun 5, 2024
8 checks passed

MFraters mentioned this pull request Jun 6, 2024

Add changelog entry for CPO elastic symmetry decomposition. #5840

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add cpo elastic tensor decomposition plugin #5116

Add cpo elastic tensor decomposition plugin #5116

MFraters commented Mar 21, 2023

bobmyhill commented Jul 8, 2023 •

edited

Loading

bobmyhill Jul 7, 2023

bobmyhill Jul 7, 2023

bobmyhill Jul 7, 2023

bobmyhill Jul 7, 2023

bobmyhill Jul 8, 2023

bobmyhill Jul 8, 2023

bobmyhill Jul 8, 2023

bobmyhill Jul 8, 2023

MFraters left a comment •

edited

Loading

MFraters Aug 3, 2023

MFraters Aug 3, 2023

gassmoeller left a comment

gassmoeller Oct 8, 2023

MFraters Mar 16, 2024

MFraters left a comment

MFraters Mar 16, 2024

MFraters Mar 16, 2024

MFraters commented May 29, 2024

gassmoeller left a comment

gassmoeller May 30, 2024

gassmoeller May 30, 2024

gassmoeller May 30, 2024

gassmoeller May 30, 2024

gassmoeller May 30, 2024

MFraters commented Jun 4, 2024

gassmoeller commented Jun 4, 2024

MFraters commented Jun 4, 2024

gassmoeller left a comment

gassmoeller commented Jun 5, 2024

MFraters commented Jun 5, 2024

	* Rotate a 3D 4th order tensor with an other 3D 4th
	* Rotate a 3x3x3x3 elastic tensor with a 3x3 tensor

	* Rotate a 6x6 voigt matrix with an other 3D 4th
	* Rotate a 4th order elastic tensor in Voigt form (6x6) with a 3x3 rotation tensor.

Add cpo elastic tensor decomposition plugin #5116

Add cpo elastic tensor decomposition plugin #5116

Conversation

MFraters commented Mar 21, 2023

bobmyhill commented Jul 8, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters commented May 29, 2024

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters commented Jun 4, 2024

gassmoeller commented Jun 4, 2024

MFraters commented Jun 4, 2024

gassmoeller left a comment

Choose a reason for hiding this comment

gassmoeller commented Jun 5, 2024

MFraters commented Jun 5, 2024

bobmyhill commented Jul 8, 2023 •

edited

Loading

MFraters left a comment •

edited

Loading