Entropy initial temperature lookup

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

In this PR, I create a new model in the initial composition to convert the initial temperature to the initial entropy by looking up a table. In the table, the pressure and the temperature need to be assigned as X and Y axis, respectively, and the entropy needs to be looked up.
After this PR, one could provide an initial temperature field but use the entropy method for computation.
This is currently working and near completion.
I want to ask for some help and review from @jdannberg, @RanpengLi,@bobmyhill before wrapping up (listed below).
FYI @mibillen

What I have done:

• A new "entropy table lookup" model in the initial composition

• A lookup table for pyrolite that does (pressure, temperature) -> entropy

• A new test entropy_initial_lookup

This test takes a constant temperature field (left subfigure) and converts it to an entropy field decreasing with depth (right subfigure).
The values of entropy fit the look-up table I provide.

Update part 1:

1. I generated a (P, S) -> T table from the (P, T) -> S table I have for pyrolite.
2. I made use of this new table in the "entropy model" Material model. So I should have two paired table and the initial temperature is converted to entropy and then back to temperature at the first step. The initial temperature is a constant value of 1600.1243 and the result is :

Update part 2:

• I Add another test that uses the WorldBuilder as the initial temperature model and generates an initial slab (same as used for the phase transition cookbook). The results for step 0 is attached:

Figure for case with an initial slab temperature field from the WorldBuilder:
left: initial temperature from the world build
middle: the entropy field
right: temperature after initial step (back projected from the entropy field by table looking up).

From the result, the initial entropy field is continuous and the back projection to temperature field is consistent with the initial temperature.

What I need help with before finishing:

• Discuss the implementation during user meetings.

• A consistent table that converts (pressure, entropy) -> temperature.

I am playing with a pyrolite composition. The model I generated from Perple_X is (temperature, pressure)-> entropy, so I need this one additional table to check the consistency of this approach (that I could get the same temperature from entropy).

• NaN values in the table

The values I have from Perple_X have NaN values. Currently, I assign them to a big value (10000.0), I'd like to know whether this is a good option.

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.

[gassmoeller]

Hi @lhy11009, the idea sounds interesting and we can discuss details during a user meeting. However, my initial question would be: Why do you need this plugin in addition to the ascii_data plugin? Or phrased differently: Could the ascii_data plugin not do already what you want to do?

[lhy11009]

@gassmoeller That's indeed an important question to think through. It would be very helpful to clarify what we want to achieve here.
If I understand it right, the ascii_data plugin will specify a composition field. In this case of application, we could use it for an entropy field. That's indeed very helpful. Moreover, I think what you mean is to have an initial entropy field (e.g. a slab) and use the ascii_data plugin, is that right?
If that what's the question contains, the answer is both yes and no. It could technically cover what I want to model but it's not going to achieve the same level of neatness this PR seems to provide.
The idea of this PR and the new plugin is more about changing the workflow. The new workflow will be using the WorldBuilder for the initial condition and have the (T, P) -> entropy table convert it into an initial entropy field.
In other words, we are trying to wrap up the entropy method as a black box for the user. Provided two (or two sets of) lookup tables, (T,P) -> S and (P, S) -> T, are given, the user would proceed the same way as if the entropy method is not used, with all the phase transitions magically show up.
Please tell me what you think.

[bobmyhill]

I like the idea, though I wonder whether it would be better to put new parameters in Material model > Entropy model (like, for example, set Initial entropy from temperature field = True) rather than Initial composition model > Entropy table lookup*.

Also, we should require input in SI units (you have currently hardcoded pressure in bar).

I'll leave code-structure comments to @gassmoeller, @jdannberg and @RanpengLi as this work will dovetail with what they are working on.

* I already find it unintuitive to have Entropy model as a Material model rather than as a different model class.

[lhy11009]

set Initial entropy from temperature field = True
This is a good idea, I could easily make the change.

we should require input in SI units

For the current version of the PR, I copied a bit from the plugin in the entropy_adiabat benchmark and used their convention. So I could make the change after a discussion with @jdannberg and @RanpengLi

I already find it unintuitive to have Entropy model as a Material model rather than as a different model class.

I agree, if this could be somehow managed, it will make the parameters I need in the PR more intuitive (don't need to be in the composition section)

I'll try to participate in the next user meeting (08/16) and discuss.

[gassmoeller]

@lhy11009 I will leave a more detailed comment later, just wanted to let you know I cannot make the user meeting tomorrow. We could either meet separately, or during the next user meeting, whichever you prefer.

[gassmoeller]

Some comments on the implementation, generally this looks good. Forget my comment about using the EntropyReader class from #5341, you need the inverse operation here, and your implementation is correct for that.

As you mentioned and as we discussed today the main challenge for you right now is that to test your implementation (and to run practical models) you need tables that allow you to go from temperature to entropy and back to temperature for the same composition. Let me know when you have figured out a good testcase then we can proceed with this PR.

[lhy11009]

@gassmoeller. I have run into an issue in testing. I have put the plugins I need into a separate .cc file that carries the name of the test, but the run fails because of it. It seems that the compilation of that .cc file is wrong without locating the right mpi.h file. I took a step back and ran the entropy_adiabat test and it produced the same error. Here is the error output when running the entropy_adiabat test (aspect/build_entropy_initial_temperature/Testing/Temporary/LastTest.log).
LastTest.log

Can you help me figure out what the issue is?

[lhy11009]

@gassmoeller , I think this PR is already in good shape. We could proceed and start the reviewing. There is one remaining issue I need help:

• run the test with the .cc file with the same name (see above).

And some task to do after it:

• reset the option to set Initial entropy from temperature field = True

Maybe we can meet in the next user meeting.

[gassmoeller]

Can you help me figure out what the issue is?

I am not sure what the error on your system is, but the one on the reference tester looks different (see https://github.com/geodynamics/aspect/actions/runs/6065701057/job/16455501549?pr=5360#step:6:607). I think you have a conflict with some changes that happened in #5365 and one of the files you include in your test is no longer a separate plugin but instead moved into main aspect. Can you rebase your PR, remove the include line that causes the error, and try again?

And we can certainly talk next week in the user meeting about next steps.

[lhy11009]

@gassmoeller You are right, there is some different between mine and the tester machine. I rebased my branch and it doesn't resolve this issue. The running of the entropy_adiabat test still crushes with the same error message. I have to switch to using the docker image for generating test results. The running of the entropy_adiabat succeeds within the docker image, demostrating the issue is indeed the difference from the tester machine.
I do have some further question about this approach:

• Is the file in the docker somehow linked to the file in my local directory. In this case, if I adjust the contents in my test prm file, is the change updated on the docker side as well? Or should I finalize my test prm file before running the docker.
• If I run the "ASPECT_GENERATE_REFERENCE_OUTPUT=1 ctest -j 4" command to generate the files in the docker, would that update my local test files?

We can wrap these questions with discussion of the further steps during our next meeting. I have a conflict of time with the next user meeting. Does it work if we schedule another meeting? Thanks for your help.

[gassmoeller]

Is the file in the docker somehow linked to the file in my local directory. In this case, if I adjust the contents in my test prm file, is the change updated on the docker side as well? Or should I finalize my test prm file before running the docker.

Files outside and inside the container are linked if you followed the instructions (i.e. if you use docker run [...] -v [...]:/opt/aspect then /opt/aspect inside the container and your aspect directory outside of the container are the same directory).

If I run the "ASPECT_GENERATE_REFERENCE_OUTPUT=1 ctest -j 4" command to generate the files in the docker, would that update my local test files?

Yes.

Based on the files you sent I cannot see why the test on your system would fail, but would work inside the container. One thing you can try is deleting the tests/ folder in the build directory, then rerun make setup_tests and rerun the test. Maybe it helps. If not we can meet sometime next week (best Wednesday or Thursday) and you can show me how you run the test.

[lhy11009]

@gassmoeller I indented the files and fix the test outputs. Now it's ready for review.

[gassmoeller]

Apologies for the long delay to review :-(. Very nice functionality, almost ready to merge. Please run make indent one more time and address the small comments I left.

[lhy11009]

@gassmoeller, I have addressed the comments and made indent.