Fix multiple phase transitions

[gassmoeller]

Hi @lhy11009 and @mibillen. I may have found a way to modify the phase averaging function to work correctly even if property values are in the .prm file and some phase transitions are only active for certain temperature ranges. The trick is in 8412483. Instead of computing the difference across the phase transition as phase_i+1 - phase_i, I instead compute it as phase i+1 - current_average. current_average is the computed property value including all previous phase transitions, but only the ones that are active for this temperature/pressure condition. I have tested this by modifying your tests for basalt and harzburgite to have values in the input parameter file and still achieve the same results as before. I also tested all existing tests for phase transitions and they still seem to work correctly. However, I have not spend a lot of time thinking about possible use cases that could go wrong with this new formulation. Can you take a look and test if my proposed change here works for your use cases? If this works then maybe geodynamics#4868 will not even be necessary.

[lhy11009]

Hi @gassmoeller . I tried to understand the change but I didn't quite figure out why it works. However, it works perfectly! Here are the plots of the 3 compositions:

Results for the basalt:

Results for the harzburgite:

Results for the pyrolite:

Maybe you can explain this to us during next week's user meeting?
Given the results, I think there are a few things I could do in the PR:

1. Merge your change with this PR (or should we put this change in another one?) It seems fine to have this in the current PR.
2. Change the tests to using this new change and use the values of densities for every phase.
3. Compose a cookbook in another PR using this new method of Phase transition.
   What do you think?

[lhy11009]

A conceptual question is, if I understand the current implementation partially, some of the phase transitions be skipped while the properties would still be the correct value related to the P, T condition? Is that right? In other words, to change for phase 0 to 4, it's not necessary that all the phase functions have to be 1.0, actually, only the last one needs to be one.

[lhy11009]

Hi Rene

I want to see if I understand the current implementation. I am using the eclogite transition as an example. So I am using two phases here and with 3 phase transitions. It seems that the new implementation is right. I draw a red path connecting the two phases. In this case, only the 2nd phase transition is in effect with the red path I draw.

Maybe we can meet during the user meeting.

[gassmoeller]

Hi Haoyuan,
Yes let's discuss this at the user meeting, but I think you already got the main idea. Btw your figure looks almost exactly as the one I drew on my blackboard :-).
The main change is that when computing the density difference across a phase transition, we no longer rely on the difference between the new phase and the old phase (including all the density changes of previous transitions), but instead we use the density difference between the new phase and the density computed so far (therefore excluding the phase transitions that didnt happen because we were outside their temperature bounds).

Glad to hear it works. Feel free to merge this PR or to just include my change into your PR.

[lhy11009]

Here is my derivation of the interaction between the two phase transitions. In this case, the second phase is not active. The phase transitions follow the red line (phase 1 -> 3 -> 4).