CF Merger #3: Compositional Flow models

[vlipovac]

PR for the next step in including the branch pmgbergen/porepy/tree/composite-flow

This PR contains the compositional flow model, related mixins, as well as a 2-component tracer flow example with tutorial.

Proposed changes

The file compositional_flow implements the new model and is the main contribution of this PR, and I think it's structure is self-explanatory as to which model parts are introduced and required for a general CF model.

The file initial_condition introduces a new layer for easier definition of IC in a multi-physics setting with multiple variables. The layer is added to existing models.

A new type of equation is added to abstract_equations, called LocalElimination. It allows to eliminate some variable using a function and some declared dependencies (other variables). The framework handles both IC and BC for the eliminated variable and the SurrogateFactory is used in combination with the function argument to obtain an AD representation of this local elimination.

The file tracer_flow contains a single-phase, 2-component flow model (tracer flow), which is showcasing how to set up a simple CF model. A new tutorial is accompanying this model and explains the steps in more detail.

The sub-package porepy/examples/geothermal_flow contains the example models provided by Omar and Michael. Due to my parental leave I am not aware how much you have discussed about how to include them. I put them in this draft since they serve as excellent examples on what is required to set up a CF model.

NOTE 28-11-2024: After discussion in person, the geothermal_flow will be removed from this PR.

Testing

@mikeljordan proposed a test using a linear tracer flow model (single-phase, 2-component) in combination with an analytical solution. The respective test can be found here here

A second tests extends the first test artificially to include additional 2 phases with fixed saturations of 1/3 and identical IC, BC, densities and viscosities.
The model per se has no meaning, but it tests the full machinery of multi-phase, multi-component flow, including having phase properties as surrogate operators and local equation framework.
Additionally, the two-variable energy balance is introduced to check whether the model is able to reflect isothermal conditions, adding an additional verification layer.

Unification of notion of fluid

The CF models have a generalized notion of fluid allowing multiple phases and components, currently coded as a model attribute fluid_mixture, based on the recently introduced extension porepy.compositional.
The (legacy) models have a data structure for fluid constants, on which various constitutive laws are built.
After talking with @keileg , the idea is for me to go ahead and propose a way to wrap the existing fluid constants into a generic 1-phase-1-component mixture, which then hopefully will be integrate into the existing models.
It is also worthwhile to consider extending the model attribute fluid to be a general FluidMixture, in order to not have two various concepts

NOTE 31-10-2024: Unification is on the way in #1244
NOTE 5-11-2024: Notion of fluid is unified and scripts were adapted correspondingly (see #1252 )

Code unification and recycling with existing fluid_mass_balance and energy_balance

There is vast space for code recycling, since the CF model currently does many things on its own, which were to a less general degree implemented in the existing fluid and energy models. My sense is that with some minor refactoring, the existing models can be included in the CF.

This point is mainly connected to the notion of fluid, and how to access relevant members such as fluid density, enthalpy etc.
But also some thoughts on the advective flux are required (f.e. a general declaration of the non-linear weight, which is in different forms implemented in the existing energy_balance and fluid_mass_balance).

If I am not mistaken, @IvarStefansson this is most relevant for you.

NOTE 31-10-2024: Unification is planned. First, minor steps were taken in #1252 (see #1252 (comment)), more thought on a general AdvectionEquation or something of that sorts is required.
NOTE 28-11-2024:

1. Unification is mostly done locally on branch. Main changes revolve around the notion of mobility, which needed to be adapted in the existing MassBalanceEquations and EnergyBalanceEquations for them to be usable in the CF setting as pressure and total energy balance equations respectively. Inconsistencies remaining in the discretization of hyperbolic/advective parts.
2. The handling of equations, variables, IC and BC is getting messy with increasing model complexity. Some abstraction layers were added for IC and BC, and the Python's super() functionality is now fully utilized in coupled-physics models with multiple mixins representing equations and variables.
3. The generalization in terms of an AdvectionEquation is not clear at this point, mostly due to unclear handling of BC for elliptic and hyperbolic discretizations (MPFA & upwinding).

NOTE 09-12-2024:
A unification of equations was done by refactoring some inner workings of existing mass and energy balance. They serve now as balance equations for total mass and energy.

Types of changes

• 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

• 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).
• pytest was run with the --run-skipped flag.

[IvarStefansson]

In the tutorial,

• please define all symbols (\Phi, \rho etc).
• Since z does not have units of kg, I find the following sentence slightly confusing "The mass corresponding to a component $i$ is represented by an overall fraction $z_i$." Is that me being stupid, or can the sentence be phrased more clearly and presicely?
• Should we refer somewhere for the phrase "Upwind-coupling."?
• In the summmary of model set-up, first bullet point, does "fluid" refer to TracerFluid? If yes, please use the class name for clarity.
• strictly speaking, BC classes provide BCs for the fluxes as well as the potential. I don't know whether the bullet point on BCs and ICs needs this nuance, though. However, I think the explanation of super should be removed from the summary. There is already a comment to a similar effect above. If you feel that one isn't exhaustive, please elaborate there rather than repeating.
• Couldn't plot_schedule be defined inside after_nonlinear_convergence? I'm scared of having mutable class attributes, even if it's hard to imagine what could go seriously wrong here.
  -The simulation example is very instructive! But I'm not sure what to think about the hypothesis. In the end, it sort of only reflects badly on the modeller - it's not like the observed results are very difficult to predict.

[vlipovac]

@IvarStefansson I have reworked the tutorial according to your review, except your point about the reference to Upwind-coupling.
I am not sure what to reference to, because Upwinding is explained to various degrees in 2 or 3 tutorials.
(At some point while writing the tutorial, I wanted to write one for Upwinding until I noticed the content in the existing tutorials. Maybe this is a nice todo for the user group.)

Regarding the tutorial design and the "bad hypothesis" (which made @keileg also "sit up" and pay attention), that was to some degree my intention. I had primarily students who are using porepy in mind (for whom impermeable fractures are not that basic). If I got the audience wrong, or if you would nevertheless prefer a more templated story-telling, let me know.

Regarding your remarks from Teams:

1. The geometry is now based on SquareDomainOrthogonalFractures, with an override for the fractures.
2. I checked out the available BC (like `BoundaryConditionsMassDirNorthSouth'), but given the setup in the tutorial, it would be a hassle to inherit from those and modify them. Also, setting BC and IC is kind of an atomic task in a model setup. There is little room for code recycling, if the setup is different.

------------- EDIT:
I modified the tracer tutorial and removed the "bad hypothesis". The sole focus of the tutorial is using the code and framework (as discussed on Teams).

[keileg]

Great work. Thanks for the effort, @vlipovac and @OmarDuran!