Add particle property plugin CPO core

[MFraters]

This is part one of adding CPO to aspect (see issue #3885). I am not entirely sure what the best approach is to get this all merged, so I tried to cut out as much as possible and focus this pull request mainly on the structure and testing the storing and loading of the CPO data in to the particle data vector.

I have kept the advection algorithm, but only the most simple derivatives algorithm. Since the D-Rex algorithm implementation is a significant amount of complicated code, I left it out of this pull request, but I did keep some of the reference to it in the code. Although the world builder initial conditions coupling is not that much code, I removed it so that I don't have to add the tests in this pull request.

Although I could already make a test for the advection schemes in this pull request, it will be a lot easier to do that once the bingham average particle property plugin has been added. So I would like to wait for this part to be merged (all the other parts are dependent on this part).

If there are strong opinions on what should be included or excluded in this pull request, let me know.

[bobmyhill]

You might find aspell useful for spell checking other big chunks of code - this is what I use once a year to catch all the new typos, and it's much less painful than reading line by line.

[MFraters]

Thanks again for the review. I addressed all your comments. Based on your review, there are just two outstanding issues in my mind:

1. Should CPO in 2D be allowed?
2. And do we want to limit the sum of the volume fractions to 1. Currently I am thinking no on both, but I don't feel to strongly about it.

[bobmyhill]

Hi @MFraters, a 2D implementation would only make sense if 2D deformation in a(n incompressible) 3D geometry model results in 2D CPO development (i.e. rotations are in-plane). Is this the case for the models you're implementing? Your simulations are likely to be CPU- and memory-intensive, so I could potentially see the benefit of a 2D (incompressible) CPO model. But it also makes for more possible cases to test, so I'm fine with limiting the implementation to 3D.

From how you've described the model, each grain rotates independently. Therefore, it should be possible to postprocess the simulation to determine the CPO development based on different volume fractions (iff CPO developed passively in the model and if one wished to do that). I'd therefore favour asserting that the volume fractions sum to one, as this will be necessary for any models where CPO has an effect on the rheology.

[MFraters]

@bobmyhill I have disabled both. It will be easier to loosen restrictions than to make it stricter. I have put the asserts in the parse parameters, which I think should be sufficient.

[bobmyhill]

Great, thanks @MFraters! Could someone else take a look at this? This is a big chunk of code, and it would be good to get one of the senseis to approve :)

[MFraters]

@gassmoeller Thanks for the review! I have addressed your comments and left a few discussions unresolved. I also found that one of the variables in the advection scheme is not needed, because it can be dealt with directly in the computation of the derivatives.

[adigitoleo]

I hope this is an appropriate place for a quick comment: The rotation matrices used for the unit test seem rather arbitrary, I wonder if you have also tested a simpler configuration e.g.

     cos(-45)    0         -sin(-45)
a =   0          1          0
     -sin(-45)   0         cos(-45)

i.e. a rotation of -45 deg about the global Y-axis. I only ask because eq. 8 in the paper produces slip invariants that are all zero in some cases, for example with the above initial rotation and a simiple shear strain rate:

0  0  1
0  0  0
1  0  0

I've written up a quick example showing how the theoretical equations can produce zero for all invariants for simple rotations at certain angles (corresponding to total rigid body rotation): http://mathb.in/71613

Depending on how it is implemented, this has the potential to crash the simulation with zero-division errors (unlikely for randomized initial orientations).

[adigitoleo]

Note that in the linked derivation, I swapped the n^\nu and l^\nu vector definitions from table 2 in the paper, because that is how the original fortran seems to have been implemented, and it also seems to be what you have done in the supplementary derivation (part 4). ~~However, I'm now thinking that you were correct and we should in fact be using the transpose of the orientation matrix to extract the slip invariants I (for orientation matrices defined using the usual alibi convention). Due to the symmetry of the orientation matrix, this still produces the same results for I and gamma in the case I considered above.~~

EDIT: After more reading, I think the confusion here arises from the fact that Kaminski 2001 and the reference code use the passive rotation convention for their input (eq. 1 Kaminski 2001). So I guess things are fine. It might be worth clarifying the rotation conventions for the input in the eventual documentation.

[MFraters]

@gassmoeller Thanks for the review. I have addressed you comments.

You currently copy all properties from its storage location into new vectors in pack/unpack_data only because you can access them more conveniently that way. The same could be done by having utility functions, e.g. get_mineral_volume_fraction(i_mineral) or get_grain_volume_fraction(i_mineral, i_grain), etc.. That way only these utility functions would need to know about the data layout. Do you think it would be better to have such accessor functions and not create new vectors for all your properties?

I would agree that that is the better option, but it would be a significant amount of extra work for making, integrating and testing the new functions. Since I already have test that this works and I don't think the gain will big enough to delay this pull request any further. It can be refactored later.

@adigitoleo Thank you for your comments. There are many ways to put together Euler angles, and Euler angles also have some issues with for example gimbal lock. That is the reaons the code is written that is internally uses rotation matrices. I agree that when we allow Euler angles for input and output, it needs to be documented well, but that will happen in a follow up pull request. Let me know if I misunderstood you comments.

[MFraters]

@gassmoeller I rebased as discussed. Let me know if you have any other comments.

@bobmyhill since your review status is still on changes requested, can I request a re-review? Thanks!

[MFraters]

I am not sure what happened, because I am sure I had work on a version in which I removed the pack and unpack functions (which I thought I put here), but I haven't been able to find it. I should have noticed that, so sorry for that. I removed them again in favor of function which directly access the data array. The only data which is stored locally now is the matrix, since it is a lot easier to to the calculations with that. But that is now only done where it is needed.

I decided to remove the Crank-Nicolson since it adds quite a bit of complexity, while also being much more expensive without having much evidence that it is much more accurate in practice. If desired it can always be added back in later.

I also change the random number generator to mt19937 as suggested and made sure to do it in a separate commit, which worked before (and after) the change.

I went through all the unresolved comments, and made sure that they are resolved and set to resolved. Let me know if there are any other comments.

[MFraters]

@gassmoeller Thanks for the review. I have addressed your comments :)

[MFraters]

Thanks for the quick review! I have addressed your comments :)

[MFraters]

Thanks again for the quick review! hmm, I am fine with adding a changelog for this if you think it is important, but this pull request just adds some core infrastructure, which is not meant to be used yet. Wouldn't it be better to wait until some actual functionality is ready?

[MFraters]

:tada: