Multiphase multicomponent buoyancy discretization

[OmarDuran]

Proposed changes

This PR introduces a mixed-dimensional hybrid upwind discretization for non-isothermal (enthalpy-based) multiphase, multicomponent flow. The implementation is based on the fractional flow formulation of the compositional equations I proposed a few years ago. It reuses the existing upwind discretization to minimize changes to the framework. Support for the standard form of the equations is not included in this PR but could be added by unfolding the definitions of the fractional flow formulation.

The following test are included

porepy/tests/models/test_fluid_property_library.py
This is a unit test for FluidBuoyancy.

porepy/tests/functional/test_buoyancy_flow.py
The integration/application test includes the following combination
Fluid type: N-phase N-components, N = {2,3} (x2)
Geometry type: Fixed- and Mixed-dimensional (x2)
Dimension: d = {2,3} (x2)
Order of conservation loss (overall mass, component mass, total fluid energy): order = {2,3,4} (x3)
Total number of test = 2x2x2x3 = 24
This test is a slow test and mark skipped.

porepy/tests/functional/test_buoyancy_flow_benchmark.py
This test runs a 3D simulation of two-phase flow in a vertical column under the effect of gravity. It is parameterized to run for three different density contrasts between the two fluids. The test validates the numerical results by comparing the final saturation profile against a reference analytical solution from Hayek, M., et al. (2009),
"An analytical solution for one-dimensional, two-phase, immiscible flow in a
porous medium," published in Advances in Water Resources.

This benchmark test provides a plotting functionality. It is controlled by an environment variable. To generate plots, run:
(Linux/macOS):
export RUN_PLOTS=1
pytest --run-skipped tests/functional/test_buoyancy_flow_benchmark.py
Test plot for low density contrast:

Test plot for moderate density contrast:

Test plot for high density contrast:

This test is a slow test and mark skipped.

As an example of the functionality, here you can see an example of mixed-dimensional segregation of two phases in 3D:
https://github.com/user-attachments/assets/69d0ce8e-4e9a-45a4-a3d2-fe1be0e77b94

Types of changes

What types of changes does this PR introduce to PorePy?
Put an x in the boxes that apply.

• Minor change (e.g., dependency bumps, broken links).
• Bugfix (non-breaking change which fixes an issue).
• New feature (non-breaking change which adds functionality).
• Breaking change (fix or feature that would cause existing functionality to not work as expected).
• Testing (contribution related to testing of existing or new functionality).
• Documentation (contribution related to adding, improving, or fixing documentation).
• Maintenance (e.g., improve logic and performance, remove obsolete code).
• Other:

Checklist

Put an x in the boxes that apply or explain briefly why the box is not relevant.

• The documentation is up-to-date.
• Static typing is included in the update.
• This PR does not duplicate existing functionality.
• The update is covered by the test suite (including tests added in the PR).
• If new skipped tests have been introduced in this PR, pytest was run with the --run-skipped flag.

[keileg]

Thanks for the PR, @OmarDuran, this is a great contribution!

I will start the review today or tomorrow, but it may take a little while to digest it all. From a very brief glance, it appears that I may have some questions or suggestions regarding how the functionality is included in the multiphysics models, but I have to read more carefully to move any further there. As mentioned previously, I am prepared to help out or take care of most of any modifications to the code; let's discuss when we get there.

I see that two tests are failing since the chemicals package is not included. Quite likely we need to include as a dependency connected to #1461 - we will discuss internally, so there is no need for you to take any action.

I can also have a look at what happens with test_example_params, I believe I have a fairly good idea of what is going on there.

[OmarDuran]

Thanks for the PR, @OmarDuran, this is a great contribution!

I will start the review today or tomorrow, but it may take a little while to digest it all. From a very brief glance, it appears that I may have some questions or suggestions regarding how the functionality is included in the multiphysics models, but I have to read more carefully to move any further there. As mentioned previously, I am prepared to help out or take care of most of any modifications to the code; let's discuss when we get there.

I see that two tests are failing since the chemicals package is not included. Quite likely we need to include as a dependency connected to #1461 - we will discuss internally, so there is no need for you to take any action.

I can also have a look at what happens with test_example_params, I believe I have a fairly good idea of what is going on there.

Regarding the test failures, they are now resolved. The dependency on chemicals is important for practical applications, so I would appreciate it if #1461 could be merged soon. For this PR, however, I have refactored and removed the dependency on chemicals.

For further modifications, I would like to stay in the loop, as I am finalizing a manuscript on this subject and need to closely follow any changes.

Regarding the inclusion in multiphysics models. The main challenge, in my view, is distinguishing what belongs to the discretization from what belongs to the reformulated expressions. I also find it difficult to understand why FluidMobility is classified as a fluid property. In the absence of a better place, I also placed FluidBuoyancy there, but in my opinion, a new module should be created fractional_flow_functions.py. and place there FluidMobility and FluidBuoyancy.

[keileg]

Thanks for taking care of the tests. We also hope #1461 to be completed soon, though it will still take a little while due to other commitments.

For further modifications, I would like to stay in the loop, as I am finalizing a manuscript on this subject and need to closely follow any changes.

Of course!

The main challenge, in my view, is distinguishing what belongs to the discretization from what belongs to the reformulated expressions.

This is indeed a good question, that part of the code is best though of as still under construction. We may revisit this as part of this PR, but it could be better to do that in a separate process to limit the number of moving parts. I will have to read the code more deeply to make an assessment.

[keileg]

I think this mainly looks good. The discretization itself is very well structured and should be easy to follow. The few questions that I have are minor stuff, partly on documentation.

I pushed a change that should speed up the update of some parameter dictionaries by avoiding Ad evaluation inside a tight loop. The computations remain really slow, I guess due to the size of the Ad operator tree for upwinding. To try to be constructive, this should be a good test case for improving evaluation speed.

Location of the upwinding: I agree that fluid_property_library is not a logical home for this code. The most natural option is constitutive_laws, but as you know all too well, that file is far too large already. I suggest we leave this for now and then possibly move later - it should make no difference to your future use of the code (right?).

To move forward, I suggest that you have a look at my comments and see what makes sense and not. Let me know if you want me to take care of the implementation of some of the parts. I would also like to add some more documentation, including method docstrings, so that I will recall what the methods are doing in half a year. I'll try to do this over the next 1-2 days, if you don't mind; no code will be touched, so there should be no impact for your work.

Let me know what you think.

[OmarDuran]

I think this mainly looks good. The discretization itself is very well structured and should be easy to follow. The few questions that I have are minor stuff, partly on documentation.

I pushed a change that should speed up the update of some parameter dictionaries by avoiding Ad evaluation inside a tight loop. The computations remain really slow, I guess due to the size of the Ad operator tree for upwinding. To try to be constructive, this should be a good test case for improving evaluation speed.

I think the following context is important. This is a non-isothermal, compositional multiphase implementation that relies heavily on the Upwind and UpwindCoupling classes. The test and benchmark files run simulations with at most three components and three phases, and as you know, nonlinear mixed-dimensional simulations in PorePy are naturally slow.

I can improve the runtime by roughly half, but the current implementation is already as fast as possible given the extra upwinding. Further speed-ups in the Upwind and UpwindCoupling classes would require clear guidance from your end on how to do so. The good news is that if a faster implementation of the Upwind and UpwindCoupling classes is available, this code will automatically benefit from it.

Another point is that including all terms required for buoyancy effects in non-isothermal compositional multiphase equations naturally increases the AD operator tree, and including capillarity would make it even larger. If you expand and account all those non-linear terms, it makes this clear.

Overall, the proposed implementation is compact and can be reused for capillary effects, which is a major advantage.

Location of the upwinding: I agree that fluid_property_library is not a logical home for this code. The most natural option is constitutive_laws, but as you know all too well, that file is far too large already. I suggest we leave this for now and then possibly move later - it should make no difference to your future use of the code (right?).

To be honest, I also don’t find it logical to mix constitutive laws with discretizations. I think we can follow up this in another PR if needed.

To move forward, I suggest that you have a look at my comments and see what makes sense and not. Let me know if you want me to take care of the implementation of some of the parts. I would also like to add some more documentation, including method docstrings, so that I will recall what the methods are doing in half a year. I'll try to do this over the next 1-2 days, if you don't mind; no code will be touched, so there should be no impact for your work.

I have almost addressed your comments.

Let me know what you think.

[keileg]

To clarify my comment regarding computational cost: This was not meant as criticism of the implementation - the formulation is what it is and as you point out, the assembly scales poorly with growing size of the operator trees. Rather, my point was that the availability of this functionality will accentuate the need for doing more fundamental work on the assembly, but that is my job, not yours. So, no worries!

[keileg]

I pushed some changes:

1. Docstrings are added to the upwinding discretization. I believe the description is fair and there is no changes to the code.
2. The reference arrays used in test_buoyancy_flow_benchmark.py are moved to applications.test_utils together with similar arrays.
3. To limit the computational cost for test_buoyancy_flow.py, I removed (commented out) the tests run for lower order of Newton accuracy and for non-fractured domains. My thinking is that the tests that are left cover everything we want from the discretization, and that it is very hard to imagine a bug that is not uncovered by the active test, but would have been detected by those commented out.
4. The benchmark test was too slow to be included even among the skipped test, so I change the name so that pytest will not detect it. It should still be used to verify the buoyancy discretization when relevant, but it is too cumbersome to use as a general skipped test.

Does this sound fair to you?

[OmarDuran]

I pushed some changes:

1. Docstrings are added to the upwinding discretization. I believe the description is fair and there is no changes to the code.
2. The reference arrays used in test_buoyancy_flow_benchmark.py are moved to applications.test_utils together with similar arrays.
3. To limit the computational cost for test_buoyancy_flow.py, I removed (commented out) the tests run for lower order of Newton accuracy and for non-fractured domains. My thinking is that the tests that are left cover everything we want from the discretization, and that it is very hard to imagine a bug that is not uncovered by the active test, but would have been detected by those commented out.
4. The benchmark test was too slow to be included even among the skipped test, so I change the name so that pytest will not detect it. It should still be used to verify the buoyancy discretization when relevant, but it is too cumbersome to use as a general skipped test.

Does this sound fair to you?

Thanks for the updates.

The docstrings for the upwinding discretization look good, no issues there. I agree with reducing the test load in test_buoyancy_flow.py; I think we should delete what is not used to avoid commented-out code. The benchmark test was a way to demonstrate the validity of the approach, and I am fine with your changes there.

Overall, the changes look good to me.

[keileg]

Thanks. I will delete the commented lines and then approve and merge.