Discrete phase function

[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.

Headline

I added a new class of phase function that handles discrete phase transitions by looking up the most dominant phases in a lookup table.

This is now ready for a review. @gassmoeller technically advises me in the PR. @ryanstoner1, I borrowed some of the ideas from your previous PR. I want to thank both of your perspectives. I also want to keep @mibillen posted on the progress of this work. This is really cool idea. To those who are interested, this PR would add a flexible utility to look up phase indicators in the material data files and use a distinct set of rheologic parameters. The new class we defined here only requires a small tweak in the material model itself, and thus would be portable to material models other than the visco plastic model.

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]

Implementation Overview

• The addition of a new phase function class, PhaseFunctionDiscrete. All the methods from the PhaseFunction class are being used in this new class, and no new methods are required. This could allow us to consider merging them with a virtual base class later on.
• Within PhaseFunctionDiscrete, a list of pointers for StructuredDataLookup is stored and loaded with a material data file upon initialization.
• Each material file name needs to be provided for every composition in the model. In the implementation, I used the variable "maph" as shorthand for "the most abundant phase" within the functions.
• The DiscretePhaseFunction class computes the value of the phase function for each phase index by determining the compositional index and referencing the corresponding lookup table. The most abundant phase is processed by assigning a value of 1.0 in the phase functions to facilitate phase transitions.
• On the other hand, the derivative of this class is set to raise error to prevent it from being called.
• I also added a function, get_number_of_coordinates, to the StructuredDataLookup class to retrieve the number of points from the T and P dimensions.
• In the ViscoPlastic class, I introduced an explicit Initialize function to load the material data file if required by the user. In its evaluate function, I simply create another instantiation of the phase function when calculating viscosity and pass it to rheology->calculate_isostrain_viscosities.

Tests

Test visco_plastic_phase_discrete

I've added a new test called visco_plastic_phase_discrete, where I use a small lookup table to conduct the test. This lookup table includes four points, with the dominant phase index assigned a value of 0 for the two low-pressure points and a value of 1 for the high-pressure points. For the rheologic parameters, I’m using diffusion creep, assigning a weaker rheology to phase 1. The results I’ve generated include pressure (P), temperature (T), and the resulting viscosity. The change in viscosity occurs around (P = 2 \times 10^{10}), which lies between the two pressure values (0 and ($4 \times 10^{10}$)) specified in the lookup table.

[lhy11009]

Test visco_plastic_phase_discrete_comp

I’ve also included a test case visco_plastic_phases_discrete_comp for the second composition, where I assigned a viscosity of ($1 \times 10^{24}$) Pa·s to the phase 0 of the second composition and ($1 \times 10^{14}$) Pa·s to the phase 1 of the same composition.

Test visco_plastic_phase_discrete_comp_invert

Furthermore, I've tested invert the phase lookup, such that we look for index 0 as the most abundant phase and assign it a smaller viscosity. This demonstrates the flexibility of the implementation: arbitrary phases can be assigned distinct rheology as long as they are marked in the look-up table in the first place.

[lhy11009]

Test visco_plastic_phase_discrete_mixed

This approach is followed by another test visco_plastic_phases_discrete_mixed, assigning a 0.5, 0.5 compositional fraction to the background and the "spharz" composition to facilitate compositional mixing. The results of the mixing are shown in the attached figure. Note that the range of values differs from the previous plot.

In the region of phase 0, the background composition has a viscosity of ($1 \times 10^{20}$) Pa·s, while the "spharz" composition has a viscosity of ($1 \times 10^{24}$) Pa·s, resulting in a mixed viscosity of around ($2 \times 10^{20}$) Pa·s. In contrast, in the region of phase 1, the background composition has a viscosity of ($1 \times 10^{18}$) Pa·s, while the "spharz" composition has a viscosity of ($1 \times 10^{14}$) Pa·s, leading to a mixed viscosity of ($2 \times 10^{14}$) Pa·s.

[lhy11009]

Some related PR ideas

• Add post-process outputs for the most abundant phases
• Incorporate the visco-plastic rheology in the entropy model using this class

[gassmoeller]

This is very cool, I am happy this could be implemented so cleanly / without affecting other parts of the material model. I have some suggestions on function names and code improvements though. Let me know when you have implemented them and I can take another look.

[lhy11009]

Hi @gassmoeller, Thank you for your comments to make it in better shape. I have now addressed these comments, so you can start another review.

[gassmoeller]

Could you check why the tests are failing?

[lhy11009]

Hi Rene (@gassmoeller ), can you take a look at the issue in "no_unity-build" test. Building on my local PC is fine. Moreover, this issue is from the reference to the functions in the newly defined PhaseFunctionDiscrete class. If it raises issue there, I am confused why the link to the previous PhaseFunction class is established in the first place.

[gassmoeller]

I'll try to get to the review tomorrow, but about the test failures: It looks like you are missing an explicit instantiation of your new template class. Check the bottom of material_models/utilities.cc and what we do in line 1543 for PhaseFunction. You can check if it works on your local machine if you rerun cmake with -DASPECT_UNITY_BUILD=OFF and then recompile aspect.

[gassmoeller]

Yes, this is much closer to being ready. I have another round of comments about documentation and some data types. I mostly want to make sure someone who is not familiar with what your new phase function does can understand how to use it from the included documentation (because even for me it is a bit hard to understand some of the details with the current documentation).

[gassmoeller]

This variable needs more explanation here and also in the declare_parameters function, because I do not understand from the documentation:

• what does the length of the vector mean (i.e. how many entries will this vector have?)
• what does each entry of the vector mean? And why are they double? I would have expected unsigned int or maybe std::string as entries. Indicators as in identifying a certain phase transition by number should be an unsigned int.

[lhy11009]

Thanks for the feedback! I wanted to clarify a bit about the indicator array and how it’s set up.

The main point here is to let us selectively assign different rheologies to specific phases, rather than having a unique rheology for each phase in the table. So, if the table has phases 0, 1, and 2, and we only want a distinct rheology for phase 2, we’d just include phase 2 in the indicator array in the .prm file.

I realize this setup may seem less intuitive since it might feel clearer if each phase had its own rheology. But going this route actually avoids a lot of fine-tuning in the table itself. For example, if we want a special rheology for the lower mantle, we’d otherwise need to merge those lower mantle phases into a single value. Here, we can mostly leave the table as-is and make any specific adjustments directly in the .prm.

This way, if we need to combine multiple phases (e.g., ten phases), we can do it easily in the .prm file. And if other phases just need the default olivine rheology, we can leave them out of the indicator array altogether. This keeps things flexible and allows for a range of uses.

As for reading the phases in as doubles—it’s just a practical choice since all the other columns in the table are doubles, even though that might seem a bit odd. When generating the table, it was tricky to keep one column as integers while everything else was in double format.

Hope this clears things up! Let me know if you have more questions.

[lhy11009]

Maybe I'll still make the index unsigned int, but I guess using a double variable in the prm file would still work if the numbers in the table are integers, it that correct?

[gassmoeller]

Theoretically you can make this work with either double or unsigned int as long as you always convert correctly. But converting is always error-prone, and unsigned int is just the more natural data type for a variable that represents an integer number (i.e. a quantity that has discrete and separate values).

I would do the following:

• use unsigned int for the dominant phase index and the data from the table and prm file everywhere
• the only place where you have to convert from double to int should then be here, because the StructuredData class unfortunately interprets all the data columns as if they are double.

[lhy11009]

Yes, I have checked this is how the implementation currently works

[gassmoeller]

Ok, good that it works with unsigned int now.

Please add the following to the documentation here: For a description of the use of the phase indicators, please see the documentation of the input parameter 'Phase transition indicators' in the function declare_parameters().

[lhy11009]

Thanks, Rene, for a lot of insightful review points. Especially about the unique_ptr, I think this time I have a feel for that. The time I changed to unique_ptr, I had to go over a few bugs already. In a way, this forces me to be careful. Here, I leave 3 open conversations for you to respond to double-check.

Let me know what you think of them.

[gassmoeller]

Is this ready for another look or do you plan to make more changes first?

[lhy11009]

Rene, it's really for another look. Sorry, I forgot to put the comments in the end.

[gassmoeller]

Sorry it took so long to get back to this PR. This looks good, I just have 4 very minor comments left. Afterwards this can be merged. Thanks for your patience, I think this is a very nice addition to the material models.

[gassmoeller]

Ok, good that it works with unsigned int now.

Please add the following to the documentation here: For a description of the use of the phase indicators, please see the documentation of the input parameter 'Phase transition indicators' in the function declare_parameters().

[gassmoeller]

/rebuild

[lhy11009]

Hi Rene, I addressed these comments and rebased them into one commit.

[gassmoeller]

Great, thank you! Good to go.