Input phase properties by value differences

[lhy11009]

Pull Request Checklist. Please read and check each box with an X. Delete any part not applicable. Ask on the forum if you need help with any step.

Describe what you did in this PR and why you did it.

Before your first pull request:

• I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[lhy11009]

In this PR, I handled the change of the inputting method of properties (e.g. density, viscosity) on phases. Now, the properties could both be assigned by values or by values of differences. @jdannberg and @anne-glerum , thanks for your advice in the discussion. Would you like to take this PR?

Description

Predating this PR, entries of the phase properties are by the exact values. On the other hand, the implementation of the phase properties is actually calculating the differences between adjacent phases. Based on these values of differences, the implementation will then calculate the values of these properties from the phase-transition functions. In this way, it makes sense to have the interface directly import the values of the differences.

Tests

Test 1: input_phase_properties_by_jumps: 4 phase transitions in 2 compositions

In this test, I use the entries of differences in values between adjacent phases, this is controlled by the option of:

set Define phase properties by differences instead of exact values = true

The result of this test shows a density change across the 6 sub-regions in the model domain. Here the surface of the domain is on top (y = 100 km) and the temperature field is set to a constant value of 273 K. The left and the right sides are two different compositions. Both of them have two-phase jumps in their sub-domain and would induce a density increase of $100 kg/m^3$ with each phase transition. This is set by:

set Densities = background:3300|100|100, right:3600|100|100

Where the values of the differences between adjacent phases are used instead of the exact values on each specific phase. Below is a visualization of the result of density.

The viscosity entries, on the other hand, are handled slightly differently. The Prefactors in the flow lows are inputted by the differences in their log10 values. For example, the meaning of "1" in the string below is that the prefactor would increase by 10-fold, and the viscosity would decrease by 10-fold across a phase transition.

set Prefactors for diffusion creep = background:5.e-23|1|1, right:5.e-18

The results of the viscosity are shown below. The composition on the left (i.e. the background) decreases in viscosity by 10-fold across each of the two phase-transition. The composition on the right, on the other hand, would remain constant, as only one value is given to the prefactors.:

[lhy11009]

This may seem confusing at the first sight. I added a function "parse_map_differences_to_values" in "material_model/utilities.cc" file to deal with the new ways of input - take the inputs of differences and compute the real values. For this, I need to pass an additional entry into the "parse_map_to_double_array" function. I also added two ways of inputting these values: a. by arithmetic or b. by logarithmic. Anyway, write to me if you have any questions. Or we can first discuss it during a user meeting before you review it. Tell me what you think.

[lhy11009]

Hi, would anyone jump in and take this review. Thanks for your help.

[anne-glerum]

Hi @lhy11009 thank you for the PR. We're preparing a new ASPECT release, so the focus lies on PR's with bug fixes at the moment. But we'll make sure to come back to this PR.

[lhy11009]

@anne-glerum. By any chance, do we have reviewers available for this PR? No hurry, I am just commenting here in case it's been forgotten.

[anne-glerum]

Hi @lhy11009 , thank you for the PR. It's in good shape. I have four overall remarks & points of discussion:

1. The name of the parameter Define phase properties by exact values instead of differences is to me less intuitive than Define phase properties by differences instead of exact values. Because, up till now, we have used the exact values.
2. The parse function can keep allow_missing_keys as last argument, so that you don't need to give it as argument each time you call the function.
3. Do all the phase parameters need the option to be passed as differences? If I remember correctly, we discussed the option for the density and there I think it makes sense. But for friction angles for example, do you see a use case?
4. Why did you choose for some parameters to add logarithmically and for some not? This needs comments. Also, should it be user-specified?

[lhy11009]

@anne-glerum, a gentle reminder, do you have time to review this again?

[anne-glerum]

@lhy11009, my apologies for the long wait.
There were still some comments open, and I've added some others. Mainly, there is a discrepancy between the name in the prm file Define phase properties by differences instead of exact values and the parameter name in the code use_values_instead_of_differences. I think the parameter name in the code should be reversed, this would then have to be changed everywhere (I have not pointed out all the locations).
I also think the 'quotient' difference should be documented better.

@jdannberg I would appreciate you having a look as well, as we three discussed this change.

[lhy11009]

Thanks again, @anne-glerum. I have provided a review for every point you addressed. Meanwhile, I added some description of this piece of contribution and of the test I generated with it. Please take a second look when you have time and get back to me.

[lhy11009]

@anne-glerum Feel free to point out more issue to address. I also added some description at the start of the PR, as well as a visualization of the test results. I also plan to open up a cookbook for the composition-dependent phase transition after this one is merged. At that point, this approach would be very clear to all users and it would be straight-forward to apply the major phase transitions in the MTZ.

[lhy11009]

Hi Anne, Sorry for the previous individual comments in the conversation. You must have received a lot of emails. Next time, I know how to use the "file change" for reviewing. Do you still find anything funky in this one?

[anne-glerum]

This is looking good, thanks for the updates! I found one tiny bug, the rest of my comments are edits of the text and its formatting.

[anne-glerum]

Suggested change

	Otherwise, the phase properties are entered by M values, each for a phase.
	Otherwise, the phase properties are entered by M values, one for each phase.

[anne-glerum]

Suggested change

	Also note that when the option is true, the differences of the log values rather than the differences of the values are used
	Also note that when the option is true, the differences of the log values rather than
	the differences of the values are used

[anne-glerum]

Suggested change

	"If it is false, the parameter file will use M values, one on each phase to"
	"If it is false, the parameter file will use M values, one for each phase to"

[anne-glerum]

Can you make each line roughly the same length and not too long? This makes it easier to read the code.

[anne-glerum]

Suggested change

	"If it is false, the parameter file will use M values, one on each phase to"
	"If it is false, the parameter file will use M values, one for each phase to"

[anne-glerum]

Please make every line roughly the same short length (like in line 184 or 186).

[anne-glerum]

I think this is a copy-paste error; shouldn't it be "Stress exponents for diffusion creep" instead of "Prefactors for diffusion creep"?

[anne-glerum]

same here as for the other rheology cc files

[anne-glerum]

same

[anne-glerum]

There is also a conflict in source/utilities.cc that needs to be resolved before this PR can be merged. You can resolve a merge conflict by rebasing to an up to date 'main' branch while on the input_phase_properties branch with git rebase main. This undoes all your changes, updates to the state of the main branch and then reapplies your changes. The reapplication will probably result in a merge conflict that you then have to fix by hand before it can continue.

[lhy11009]

Hi Anne, I have addressed your comments. But unfortunately, the "parse_map_to_double_array" function has been completely rewritten in previous PRs. So I have to change my implementation as well. Can you review this part again? I have marked the points in the code.

[lhy11009]

I found the function 'parse_map_to_double_array" has been rewritten in previous PRs. So I have to the add this attribute to the class they implemented.

[gassmoeller]

I would much prefer to have this information outside of the parse function. The reason is the following:
This function can be used to parse arbitrary maps for any input parameter in ASPECT. We tried very hard to make it this generic. However, this new property is extremely specific and makes not much sense for other properties except the ones in material models. I understand that you need a way to modify the parsing if differences rather than values are expected, but I think this can be achieved in a more natural way (see comments below).

[lhy11009]

So instead of handling the values of inputs, I can only pass this to the instance of this class. This class would then go through the interface of another function and reach the point where the inputs are actually handled (you can follow it in the code).

[gassmoeller]

Hi Haoyuan,
Apologies for rewriting the parse_map function in #5055 so you had to redo some of your work, I didnt realize this PR involved changes to that function.

I had a brief look at this PR, and it generally looks in a good shape, but I want to raise the following question first:

Is there an application where you absolutely need this new feature? Can you describe it?

The reason I would like this clarified is that in your description of this PR you say: "... therefore it makes sense to have the interface directly import the values of the differences."
However, this is not exactly what you are doing here. Instead of always importing the differences, you add an additional way to input the values (in fact two, counting the log values as well), because we cannot remove the old way of providing values. As you noticed this adds a lot of complexity to the parsing function, and the fact that the new input is optional also adds an additional input parameter to all places that want to use it. In order to justify this added complexity we would need a reason to add it.

Even if there is an reason: Could there be a structurally easier way to achieve the same result? E.g. if it is more natural for the computation to use differences rather than values, but all of our input is meant for values, could we just compute the differences internally after parsing the values? If this is impossible with the current interface of the function, could we change the interface or make a new interface for the function? One option I have been thinking about during writing #5055 was to have a function return the actual multimap instead of a flattened vector of values. This would allow the calling function to easily compute the differences (or log values) for each phase internally by themselves, without having the user provide them. This way we could also save the additional input parameter, every material model or equation of state could decide for itself if it wants to deal only with values, or only with differences, or a mix of the two for different parameters (hard-coded, instead of chosen by the user).

I apologize that I bring this up now after you have put in all the work to polish this PR. This part of the input parsing has simply grown extremely complex, and I am hesitant to add further complexity to it. Let me know your thoughts, I am happy to help with changing parse_map_to_double_array if you need additional options (if we decide to move forward in this direction with this PR, I have an idea for how to change the interface to at least not add another option to the function).

[gassmoeller]

This function is missing names for the input arguments.

[gassmoeller]

I would much prefer to have this information outside of the parse function. The reason is the following:
This function can be used to parse arbitrary maps for any input parameter in ASPECT. We tried very hard to make it this generic. However, this new property is extremely specific and makes not much sense for other properties except the ones in material models. I understand that you need a way to modify the parsing if differences rather than values are expected, but I think this can be achieved in a more natural way (see comments below).

[gassmoeller]

What is the reason for handling the prefactors differently? This needs to be documented in every material model input property that uses the log values rather than the actual values. Is there a functional reason this has to be the log value, or could you use for example the CitcomS notation, where every value is a multiplier to the initial value (i.e. prefactor = base_prefactor * multiplier_for_this_phase). This would seem easier to understand (e.g. I was at first confused that a value of 1 in the input file means change by one order of magnitude, and 0 would mean change nothing).

[gassmoeller]

Please leave a note here why this may be needed or a good idea. From the conversation in this PR I gather you have discussed that at some point with @jdannberg or @anne-glerum, but I (and other ASPECT users) were not part of that conversation, and the code itself does not seem to mention it. Without this explanation this looks like a whole lot of complication.

[lhy11009]

@gassmoeller. Thanks for your comments on the PR. We have discussed the reason to have this implementation is to have it handle a more complex scenario in the phase transitions. The current implementation we have could only handle layered phase transitions. The values entered through that interface has to be layered as the change has to go from phase 1 to 2, then to 3, in the exact order. If, say, phase 1 is transferred to phase 3 directly, then it cannot be achieved with the interface we have. If, on the other hand, we could pass in the value of differences instead, there is a way to do that, the transition from phase 1 to 3 is just another phase transition, where we can assign a clayperon slope and pass in the value differences across this phase transition.

[lhy11009]

The idea of this PR is to only change the interface of the parse_map_to_double_array function so that the value exported from it remains the same no matter if the differences or the real values are entered. And originally I thought it wouldn't be a complexity when I was just considering adding up the differences to the base phase. But it turns out every interface on the way has to be changed, which I didn't expect honestly.

[lhy11009]

@gassmoeller Anyway, maybe we can involve Magali in the conversation, as she might be able to explain things better than I do, and the three of us could have a short discussion about it? I could send an email around if you say yes. Thanks for your time.

[lhy11009]

@gassmoeller. Hi Rene. I was thinking it would not be easy to explain everything here as well as the idea you have got for the implementation, so maybe a short discussion would be the best thing to do. I would vote for a 30 mins discussion. Anyway, get back to me when you have time. Thanks a lot.

[gassmoeller]

Sure, I am happy to meet. Would you have time during the ASPECT user meeting today? Otherwise we can also meet later, can you set up a when2meet or similar to find a time? I think we can probably clarify the questions in a short time, we should just take care to document the result of our discussion in this PR so that others can see the justification for whatever conclusion we come to.

[lhy11009]

@gassmoeller @mibillen
An example of the basalt phase transition on top of this other PR: #4591
Left: density, right: viscosity.

In this example, I use 3 phase transitions for the shallower basalt -> eclogite phase transition (yellow -> purple). There are assigned limits on temperature (this other PR), and values inputted by the differences between phases.

Here is the prm file for this test:
tests_visco_plastic_phases_basalt.prm.txt

To run this test, you'll need to use my branch:
https://github.com/lhy11009/aspect/tree/jump_and_limits

@mibillen also has an idea of how to implement this change otherwise, maybe we could arrange a short meeting on the next user meeting?

[gassmoeller]

Thanks for providing the example, and yes I am happy to discuss at the next user meeting.

However, I have a question about your example: To me it looks like the same behavior should have been possible with the old way of providing densities (as values instead of differences) and the following list of densities:

set Densities = background:3000.0|3540.0|3540.0|3540.0|3613.0|3871.7

Please provide an example for why the new way of specifying properties is necessary, i.e. something that would not work with our current input options..

[gassmoeller]

As an additional question: It looks like the function parse_map_differences_to_values exactly reverts the changes of your new input options. This seems unnecessary (to first change our inputs to use differences, and then revert to values internally anyway). Did I misunderstand the purpose of this function?

[lhy11009]

@gassmoeller , maybe you click the link to the other PR? There the file is different. But the prm I append in the conversation has:
set Densities = background:3000.0|540.0|540.0|540.0|73.0|258.7

Secondly, you mentioned this option in the parse function. You are right about it. My consideration with this PR is to handle everything in the parse function, which would revert the differences to the real values with the new inputs. I was thinking about maintaining a minimum change to the other parts of the code. But now that you mentioned the interfaces of the parse function, we can handle it otherwise.

[lhy11009]

@gassmoeller , maybe you click the link to the other PR? There the file is different. But the prm I append in the conversation has:
set Densities = background:3000.0|540.0|540.0|540.0|73.0|258.7

Secondly, you mentioned this option in the parse function. You are right about it. My consideration with this PR is to handle everything in the parse function, which would revert the differences to the real values with the new inputs. I was thinking about maintaining a minimum change to the other parts of the code. But now that you mentioned the interfaces of the parse function, we can handle it otherwise.

[gassmoeller]

@gassmoeller , maybe you click the link to the other PR? There the file is different. But the prm I append in the conversation has:
set Densities = background:3000.0|540.0|540.0|540.0|73.0|258.7

No, you misunderstand me. I have seen that the list of densities you give is the list of densities in your example model. My question is: Why do we need to specify it like this (as differences)? Why would the old way of specifying the densities (as values) not have worked? If you put in the following densities and use the old way of using densities, would you not get the same result?

set Densities = background:3000.0|3540.0|3540.0|3540.0|3613.0|3871.7

[lhy11009]

@gassmoeller That's a good question. I guess it was one of our initial considerations. I have tried a couple of times to explain why it doesn't work and it's not so intuitive, it's better to show this with the example you inputted there. Let me try that first.

[lhy11009]

So in this case, if we input the value on these phases, only the first PT works (then one that shows a slope between the blue and green region around 700 in the y axis. The problem is we want to have 3 boundaries for this basalt -> eclogite transition. So in general, by inputting the value, we can only take care of layered phase transitions.

[gassmoeller]

So in this case, if we input the value on these phases, only the first PT works (then one that shows a slope between the blue and green region around 700 in the y axis. The problem is we want to have 3 boundaries for this basalt -> eclogite transition. So in general, by inputting the value, we can only take care of layered phase transitions.

But this is just a description, not an explanation. I can see from your figure that inputting the values leads to a different result, but I do not understand why. First, I do not understand why the other two branches of the first phase transition (left and right of the density anomaly around 150-200 km on the x axis) now do not exist. Second, I do not understand how inputting the full values, then computing the difference, then applying the phase function is different from inputting the difference, then applying the phase function.

As far as I understand you are saying the following two operations lead to different results:

1. Input values:

density_values = rho1,rho2,rho3
rho_combined = rho1 + phase_function2 * (rho2-rho1) + phase_function3 * (rho3-rho1)

2. Input differences:

density_differences = rho1, delta_rho2, delta_rho3
rho_combined = rho1 + phase_function2 * delta_rho2 + phase_function3 * delta_rho3

Did I misunderstand how our phase transitions work?

[lhy11009]

Yes, I think it's a good point that we can write it out with the different expressions using the two methods.
Input values:
density_values = rho1,rho2,rho3
rho_combined = rho1 + phase_function2 * (rho2-rho1) + phase_function3 * (rho3-rho1)
In this case, rho1, rho2, and rho3 are three different phases, thus they have different values. So a related question is, is a phase transition always define as a new phase. In the case of layered phase transitions, we will always have a new phase with a phase transition, that's where this input-by-value works perfectly. But if the phase transition is not layered, rho1 -> rho2 and rho2 -> rho3 could be defining the multiple phase transitions between two phases. In this case, it's the value (rho2 - rho1) and (rho3 - rho2) that has to be the same. That's the example shown in the eclogite transition.
If we input the values for the eclogite transition in the example that doesn't work, the differences in the values would disappear for some of the phase transitions we want to implement.
In terms of the equations, one way it could work for inputting the differences is
example a:
density_differences = rho1, delta_rho2, delta_rho3
and that:
delta_rho2 = delta_rho3
This is possible, as the density differences across the two phases transition would be the same.
But one way it doesn't work is for inputting the values:
density_values = rho1,rho2,rho3
and that
rho2 = rho3
This phase transition would disappear if we visualized in a phase diagram I showed. That's because it's still the differences getting computed, and in this case, this defines no change in density.

[lhy11009]

Hi Rene and Magali

I am not sure I get the right time for the user meeting. I was in the zoom link but no one showed up. Would we choose another time for discussion?

[mibillen]

Haoyuan - I was in transit during the Aspect meeting this week. Rene - did Haoyuan’s example calculation help to illustrate the issue that he is trying to resolve with the eclogite transition not being a simple linear Clapeyron slope, dP/dT, transition like olivine/pyroxene… the phase boundary has 3 different segments depending on temperature (as actually some other phase transitions at the 660 do as well)? If not, then I think it would be helpful to meet briefly via zoom. Haoyuan is proposing a change that will work for the special case of the eclogite transitions AND will also continue to work for the layered case. I also suggested to him a more minor/backward compatible change to the parsing function to make this work. I can meet next week any morning between 8-9 am if that would useful. Magali

…

On Mar 15, 2023, at 9:46 AM, lhy11009 ***@***.***> wrote: Hi Rene and Magali I am not sure I get the right time for the user meeting. I was in the zoom link but no one showed up. Would we choose another time for discussion? — Reply to this email directly, view it on GitHub <#4868 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACS64UTP2A3HFN4GIQU57ILW4HW6PANCNFSM52KNIYDQ>. You are receiving this because you were mentioned.

____________________ Professor of Geophysics Chair, EPS Graduate Program Geophysics Minor Advisor Department of Earth and Planetary Sciences Univ. of California, Davis Davis, CA 95616 2129 Earth & Physical Sciences Bldg. Office Phone: (530) 752-4169 https://magalibillen.faculty.ucdavis.edu/ Students - need to meet with me? Check available appointment times: https://calendar.app.google/49asftepWVS1AFvy9 Or e-mail me list of days/times you are available & topic of discussion. Pronouns: She/her/hers First Generation College Student: https://firstgen.ucdavis.edu/

[gassmoeller]

Haoyuan you were late to the user meeting. California and Florida change to daylight savings time at the same date, so the meeting happened at the same time as usual, but we had already finished when you joined (I only saw the email later).

Magali thanks for offering, yes I think it would be good if we could meet this week. Would Thursday 8-9 am Pacific work for both of you? (I can also do Friday if Thursday doesnt work). Let me know.

In preparation for this meeting, here are my two major concerns we need to address before we can move forward with this PR:

1. All the explanation that we settle on has to be included in the code of this PR, not just this conversation. It looks like Haoyuan discussed all this with Anne and Juliane before, and I dimly remember like I have discussed it with him in the past as well, but this PR does not include this explanation and justification, which is why we keep coming back to the same discussions. It is important to include the justification in the code itself (or test or cookbooks), because other maintainers will have to understand this in the future as well. The current code only mentions that we can do this now, but not why or how.
2. I appreciate your explanations about the not layered phase transitions and it explains why we cannot handle it as one phase transition, but I do not yet see, why we need to input differences instead of values. In particular, I am confused that in all places that this PR changes, right after reading the differences, we immediately convert them back to values using the function parse_map_differences_to_values. See for example this line. If it is possible to immediately convert one way of defining the phase transition properties into the other without actual temperature/pressure conditions, why is it not possible to input the correct values from the beginning?

[lhy11009]

@gassmoeller, the time works for me. Sorry for being late at the meeting, last week was a mess. I was late for two meetings and early for one another depending on where others are. I guess I can better explain your question with a board (or a notebook on PC).

[gassmoeller]

Sorry for being late at the meeting

No worries, time change week is always confusing.

the time works for me.

Ok, let's aim for that if it works for Magali as well.

[gassmoeller]

Please also add a test that changes the plastic parameter angle of friction and cohesion.

[lhy11009]

Some additional comments:

1. We talked about inputting the prefactors by the ratio between adjacent phases instead of by the logarithmic differences. However, to make this an option, it has to be fixed in this function:
   read_or_check_map_structure
   This is from the implementation of
   if (n_values == 1)
{
const double field_value = map.find(field_name)->second;
for (unsigned int i=1; i<n_expected_values; ++i)
{
if (options.use_differences_instead_of_values)
{
map.emplace(field_name, 0.0);
}
else
{
map.emplace(field_name, field_value);
}
}
}
   If the values are the ratio between two phases, then the default value to input is 1.0 instead of 0.0. However, this function doesn't know about that information (whether a logarithm method is needed later to compute the real value).

[lhy11009]

2. The test with the plastic properties visco_plastic_phases_viscosity_jump_plastic (The changes are in Angles of internal friction)
   This test is compiled in a similar way to the first test, with variations in the angles of internal friction
   The viscosity field for this new test:
   

[lhy11009]

3. For the additional cookbook we talked about, I am going to use the implementation in this PR as well:
   Add upper and lower limits of temperature in phase transitions #4591
   So let's plan on finishing that PR first.

[gassmoeller]

@lhy11009 do you agree that this PR can be closed now that #4591 is merged?

[lhy11009]

@gassmoeller Definitely. Thanks for the reminder.