Skip to content

Reference states#1696

Open
jwboth wants to merge 27 commits into
developfrom
reference-states
Open

Reference states#1696
jwboth wants to merge 27 commits into
developfrom
reference-states

Conversation

@jwboth

@jwboth jwboth commented Jun 19, 2026

Copy link
Copy Markdown
Contributor

Proposed changes

This PR introduces an additional state of Variable and TimeDependentDenseArray following the style of previous_iteration and previous_time_step. These now also are able to hold reference values which can be obtained in AD operator form through call of the method op.reference() for any AD operator op, e.g., model.pressure(subdomains).reference(). In addition, perturbations to reference states can be conventiently fetched by op.perturbation_from_reference(). If not set, reference values are 0. No optimization in terms of emtpy operators has been implemented but could be foreseen in near future.

The API for setting values has been extended aiming at canonical extension of current workflows (EquationSystem.set/get_variable_values and pp.set/get_solution_values). An additional keyword reference needs to be set true.

As an example (hopefully of actual value) is added utilizing perturbations from reference states: Porosity with mere rock compressibility simplifying standard two-way coupled poromechanics. It currently mainly serves for testing the potential integration of reference states in models.

NOTE: This PR does not integrate reference states in models and variables like currently proposed in #1686 (initialization). The aim is to base further work on initialization on the current PR.

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.

@IvarStefansson IvarStefansson left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Main comment: Think about precendence for reference - time step - iteration

Comment thread src/porepy/compositional/materials.py
Comment thread src/porepy/numerics/ad/ad_utils.py Outdated
"""
try:
value = data[pp.REFERENCE_SOLUTIONS][name].copy()
except KeyError as err:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are we sure we want to have fallback?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would say no, I prefer to have the code crash on me. @jwboth what do you think?

@keileg keileg Jun 19, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update after a chat with @jwboth: A default value will make it possible to take the perturbations of general operators, without having to explicitly set a zero reference on all leaves (in the operator tree sense). This seems like a nice feature to have and will be useful for the initialization. The original implementation will not work, since under parsing this may lead to arithmetic operations with incompatible operands. The solution we ended up with for now is implemented below; this is not elegant and could have been less computationally expensive, but it should do the job for now.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Which types of leaves will this be invoked for? I'm not certain I grasp all possible consequences of this.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jwboth could you please comment on this one.

Comment thread src/porepy/numerics/ad/ad_utils.py Outdated
A copy of the values stored at the passed index.

"""
if reference:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this precedence makes sense. But we need to document the rules for previous time steps and iterates and references somewhere.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Moved all get and set methods to a new module, will document at the module level there.

Comment thread src/porepy/numerics/ad/operators.py Outdated
Comment thread src/porepy/numerics/ad/operators.py Outdated
Comment thread src/porepy/numerics/ad/operators.py Outdated
Comment thread src/porepy/numerics/ad/operators.py Outdated
Comment thread src/porepy/numerics/ad/operators.py
Comment thread src/porepy/numerics/ad/operators.py Outdated
Comment thread src/porepy/numerics/ad/operators.py Outdated
@keileg keileg force-pushed the reference-states branch from f854035 to 701148d Compare June 26, 2026 11:40
@jwboth jwboth mentioned this pull request Jun 27, 2026
13 tasks
@keileg keileg force-pushed the reference-states branch from 709f05b to 2a933b4 Compare June 29, 2026 07:42
- ReferenceOperator: The operator represents a reference value and will evaluate this
reference.

When combining the mixins, their behavior is prioritized as follows:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This list, which hopefully is complete, defines the prioritization logic of various operator states. Please review carefully and keep in mind when assessing the rest of the code.

return Scalar(self.solid.porosity, "porosity")


class CompressibleRockPorosity(pp.PorePyModel):

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@IvarStefansson I have not taken a stand on what to do with this class. My initial thought is, do we really need it, or is it 'nice to have'? If its main purpose is as a testing ground for initialization routines, there is information in that.

@@ -0,0 +1,353 @@
"""This module contains helper functions for setting, getting and shifting values in the

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This module is the new home for get and set methods (moved from ad_utils). I think that making a move was the right decision, whether this is the right placement is less clear.

@keileg keileg force-pushed the reference-states branch from 94e5d80 to b06a169 Compare June 29, 2026 08:31
@keileg

keileg commented Jun 29, 2026

Copy link
Copy Markdown
Contributor

I have done a partial re-implementation of this PR. Main points:

  1. Iterate, TimeDependent and Reference operators are now mixin classes. They have been moved to a separate module, operator_states.
  2. The test suite for these mixin classes has been reworked. The test coverage should be better, as should test documentation.
  3. The get- and set-methods for variable states were moved out of ad_utils to a separate file.
  4. The rules of precedence for taking reference, previous iterate and previous time step have been revised and documented at the module level in operator states. I believe the rules are implemented consistently throughout the code base, but this should be carefully examined in the review.

@keileg keileg force-pushed the reference-states branch from fcc1843 to 597109a Compare July 2, 2026 07:08
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants