Discrete phase funcion

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

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

[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

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

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