aspect
aspect copied to clipboard
geodynamics
Input phase properties by value differences
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:
- [x] I have read the guidelines in our CONTRIBUTING.md document.
For all pull requests:
- [x] I have followed the instructions for indenting my code.
For new features/models or changes of existing features:
- [x] I have tested my new feature locally to ensure it is correct.
- [x] I have created a testcase for the new feature/benchmark in the tests/ directory.
- [x] I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.
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.:
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.
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.
@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.
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.
@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.
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.
@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.
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.
@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.
@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.
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.
@gassmoeller @mibillen An example of the basalt phase transition on top of this other PR: https://github.com/geodynamics/aspect/pull/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?
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..
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?
@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 , 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 , 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
@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.
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.
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:
- Input values:
density_values = rho1,rho2,rho3
rho_combined = rho1 + phase_function2 * (rho2-rho1) + phase_function3 * (rho3-rho1)
- 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?
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.
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?
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 https://github.com/geodynamics/aspect/pull/4868#issuecomment-1470386924, 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/
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:
- 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.
- 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?
@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).
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.
Some additional comments:
- 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_structureThis is from the implementation ofif (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).