1178 robin conditions in models

[IngridKJ]

Proposed changes

Addresses issue #1178: Allow for usage of Robin type boundary conditions with the model class setups in PorePy.

This PR includes adding a robin operator parameter to the method _combine_boundary_operators(). In addition to this, each call to _combine_boundary_operators() is refactored into new methods specific for each physical problem, distinguished by their flux name:

• combine_boundary_operators_darcy_flux()
• combine_boundary_operators_fourier_flux()
• combine_boundary_operators_mechanical_stress()

Additionally: Fixes a problem of singular matrices related to nodes with only one subcell gradient but two or more boundary faces (for Robin conditions. This was already in place for Neumann).

Practicalities

The Robin conditions are implemented in the following form in PorePy: flux + alpha * potential = G. What is needed by the user to utilize such boundary conditions is to: In bc_type-methods:

• Use the flag "rob" when setting boundary condition type, analogous to setting Dirichlet ("dir") and Neuman ("neu").
• Assign the Robin weight, alpha. The Robin weight is an attribute of the boundary condition class (.robin_weight). See the __init__-methods of BoundaryCondition/BoundaryConditionVectorial (porepy/params/bc.py) for an example of what the weight should look like.

In bc_values-methods:

• Assign the right hand side G. G has units of flux, and is set by bcvalues"flux_name"() - same as the Neumann values. As both BC type values are set by the same method, it is important to make sure the Neumann and Robin values are assigned to the correct boundary sides.

Why refactoring the combine operators call? Tailoring of each operator is now easily achieved by just overriding the wrapper method (combine_boundaryoperators"flux_name"). Note: The default Robin operator is the same as the Neumann operator.

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).
• [x] New feature (non-breaking change which adds functionality).
• [x] Breaking change (fix or feature that would cause existing functionality to not work as expected).
• [x] 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.

• [x] The documentation is up-to-date.
• [x] Static typing is included in the update.
• [x] 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.

[OmarDuran]

The question is general in nature. I believe that there is room for a unification of types and data, such as

Types:

• Dirichlet_type
• Neumann_type
• Robin_type

Data:

• Dirichlet_data
• Neumann_data
• Robin_data

[IvarStefansson]

Thank you for your contribution.

It is important to consider the following:

• Perform a test to ensure that boundary segments labeled with Robin type do not intersect with other types;
• Include a test for scalar and vector problem, for which Robin approximates the functionality of Dirichlet and Neumann;
• Include a tutorial that illustrates how the boundary condition is intended to be used in the models.

To your first point, @OmarDuran, isn't that more a test for the BoundaryCondition class, rather than this Model integration?

@IvarStefansson It depends on how you define integration. The tutorial and test for boundary segment intersection can be seen as separate PRs.

[IngridKJ]

Thank you for the quick review, @OmarDuran. I have now addressed the suggested changes and comments, and all are covered in one way or another. The latest test I wrote (checking Robin limit case of alpha = 0 equals Neumann) is perhaps a bit brute force and could need some more fancy parametrization of sorts. I am however unsure how to adapt it...

[IngridKJ]

@IvarStefansson: PR is ready for (hopefully final) review now. I addressed your comment about nonzero boundary values, and after the meeting with you and Eirik today, I have also refactored the eliminate-asymmetric-stress-tensor functionality for Neumann/Robin.