Introduce source and target spaces for operator domains

CHANGELOG.md

@@ -17,13 +17,15 @@ PR 1772: SPEED: Faster construction of grid topologies on 3d simplex grids.
 PR 1770: Simplify PR template by deferring to CONTRIBUTING.md for details on code style and conventions.
 PR 1767: Bugfix in THM manufactured setup.
 PR 1766: Fix bug in generation of Cartesian grids not anchored in the origin.
+PR 1747: Ad Operators domain, range and dof info is specified by OperatorSpace objects.
 PR 1732: Added support for evaluating restricted variable subsystems without evaluating
     the full Jacobian and then slicing columns.
 
-
-
-
 ### Breaking changes
+PR 1747: DenseAdArray, SparseAdArray must be initialized with OperatorSpace objects.
+Other Ad Operators are by default treated as mapping from and to the same grid, with one
+degree of freedom per cell. To change this default, explicit OperatorSpaces must be
+assigned.
 
 Changes to the front end (user-facing code: multiphysics models, solvers, grids, and
 similar) that require users to update their own code when upgrading.

src/porepy/applications/discretizations/flux_discretization.py

@@ -17,6 +17,10 @@ class FluxDiscretization:
     fourier_keyword: str
     """Keyword for the Fourier (energy) problem."""
 
+    nd: int
+    """Ambient dimension of the problem. Normally set by a mixin instance of
+    :class:`~porepy.models.geometry.ModelGeometry`."""
+
     def darcy_flux_discretization(
         self, subdomains: list[pp.Grid]
     ) -> pp.ad.MpfaAd | pp.ad.TpfaAd:
@@ -31,9 +35,9 @@ def darcy_flux_discretization(
         """
         scheme = self.params.get("darcy_flux_discretization", "mpfa")
         if scheme.lower() == "mpfa":
-            return pp.ad.MpfaAd(self.darcy_keyword, subdomains)
+            return pp.ad.MpfaAd(self.darcy_keyword, subdomains, nd=self.nd)
         elif scheme.lower() == "tpfa":
-            return pp.ad.TpfaAd(self.darcy_keyword, subdomains)
+            return pp.ad.TpfaAd(self.darcy_keyword, subdomains, nd=self.nd)
         else:
             msg = f"{scheme} is not a valid Darcy flux discretization scheme. "
             msg += "Use either 'tpfa' or 'mpfa'."
@@ -53,9 +57,9 @@ def fourier_flux_discretization(
         """
         scheme = self.params.get("fourier_flux_discretization", "mpfa")
         if scheme.lower() == "mpfa":
-            return pp.ad.MpfaAd(self.fourier_keyword, subdomains)
+            return pp.ad.MpfaAd(self.fourier_keyword, subdomains, nd=self.nd)
         elif scheme.lower() == "tpfa":
-            return pp.ad.TpfaAd(self.fourier_keyword, subdomains)
+            return pp.ad.TpfaAd(self.fourier_keyword, subdomains, nd=self.nd)
         else:
             msg = f"{scheme} is not a valid Fourier flux discretization scheme. "
             msg += "Use either 'tpfa' or 'mpfa'."

src/porepy/applications/test_utils/benchmarks.py

@@ -84,7 +84,10 @@ def effective_tangential_permeability(
             ].values[0][0]
             values.append(val_loc)
         return pp.wrap_as_dense_ad_array(
-            np.hstack(values), size, "effective_tangential_permeability"
+            np.hstack(values),
+            size=size,
+            name="effective_tangential_permeability",
+            grids=subdomains,
         )
 
     def effective_normal_permeability(

src/porepy/applications/test_utils/well_models.py

@@ -231,5 +231,7 @@ def well_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
         """
         size = sum(sd.num_cells for sd in subdomains)
-        permeability = pp.wrap_as_dense_ad_array(1, size, name="well permeability")
+        permeability = pp.wrap_as_dense_ad_array(
+            1, size=size, name="well permeability", grids=subdomains
+        )
         return self.isotropic_second_order_tensor(subdomains, permeability)

src/porepy/compositional/compositional_mixins.py

@@ -658,7 +658,7 @@ def overall_fraction(
         if self.fluid.num_components == 1:
 
             def fraction(domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
-                return pp.ad.Scalar(1.0, "single_feed_fraction")
+                return pp.ad.Scalar(1.0, "single_feed_fraction", domains=domains)
 
         # NOTE: If the reference component fraction is independent, below elif-clause
         # will be executed, instead of the next one.
@@ -737,7 +737,7 @@ def saturation(self, phase: Phase) -> DomainFunctionType:
         if self.fluid.num_phases == 1:
 
             def saturation(subdomains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
-                return pp.ad.Scalar(1.0, "single_phase_saturation")
+                return pp.ad.Scalar(1.0, "single_phase_saturation", domains=subdomains)
 
         # NOTE: If the reference phase is independent, below elif-clause will be
         # executed.
@@ -773,7 +773,7 @@ def phase_fraction(self, phase: Phase) -> DomainFunctionType:
         if self.fluid.num_phases == 1:
 
             def fraction(domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
-                return pp.ad.Scalar(1.0, "single_phase_fraction")
+                return pp.ad.Scalar(1.0, "single_phase_fraction", domains=domains)
 
         elif get_local_equilibrium_condition(self) is None:
 
@@ -845,7 +845,9 @@ def extended_fraction(
 
             def fraction(domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
                 return pp.ad.Scalar(
-                    1.0, "single_component_single_phase_extended_fraction"
+                    1.0,
+                    "single_component_single_phase_extended_fraction",
+                    domains=domains,
                 )
 
         # If no unified equilibrium, calling the extended fractions will raise an error.
@@ -901,7 +903,9 @@ def partial_fraction(
 
             def fraction(domains: pp.SubdomainsOrBoundaries) -> pp.ad.Operator:
                 return pp.ad.Scalar(
-                    1.0, f"single_partial_fraction_{component.name}_{phase.name}"
+                    1.0,
+                    f"single_partial_fraction_{component.name}_{phase.name}",
+                    domains=domains,
                 )
 
         # Case only 1 phase in entire mixture, partial fractions are equal to overall

src/porepy/examples/flow_benchmark_2d_case_1.py

@@ -149,7 +149,10 @@ def fracture_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """
         size = sum([sd.num_cells for sd in subdomains])
         permeability = pp.wrap_as_dense_ad_array(
-            self.solid.fracture_permeability, size, name="fracture permeability"
+            self.solid.fracture_permeability,
+            size=size,
+            name="fracture permeability",
+            grids=subdomains,
         )
         return self.isotropic_second_order_tensor(subdomains, permeability)
 

src/porepy/examples/flow_benchmark_2d_case_3.py

@@ -82,7 +82,7 @@ def fracture_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
 
         """
         if len(subdomains) == 0:
-            return pp.wrap_as_dense_ad_array(1, size=0)
+            return pp.wrap_as_dense_ad_array(1, size=0, grids=[])
         permeabilities = []
         for sd in subdomains:
             permeabilities.append(
@@ -91,7 +91,9 @@ def fracture_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
                     "m^2",
                 )
             )
-        permeability = pp.wrap_as_dense_ad_array(np.concatenate(permeabilities))
+        permeability = pp.wrap_as_dense_ad_array(
+            np.concatenate(permeabilities), grids=subdomains
+        )
         return self.isotropic_second_order_tensor(subdomains, permeability)
 
     def normal_permeability(self, interfaces: list[pp.MortarGrid]) -> pp.ad.Operator:
@@ -105,7 +107,7 @@ def normal_permeability(self, interfaces: list[pp.MortarGrid]) -> pp.ad.Operator
 
         """
         if len(interfaces) == 0:
-            return pp.wrap_as_dense_ad_array(1, size=0)
+            return pp.wrap_as_dense_ad_array(1, size=0, grids=[])
         permeabilities = []
         for intf in interfaces:
             # Get hold of the fracture subdomain.
@@ -133,7 +135,7 @@ def normal_permeability(self, interfaces: list[pp.MortarGrid]) -> pp.ad.Operator
                 self.units.convert_units(val * np.ones(intf.num_cells), "m^2")
             )
         return pp.wrap_as_dense_ad_array(
-            np.hstack(permeabilities), name="normal_permeability"
+            np.hstack(permeabilities), name="normal_permeability", grids=interfaces
         )
 
 

src/porepy/examples/flow_benchmark_2d_case_4.py

@@ -142,7 +142,10 @@ def fracture_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
         """
         size = sum([sd.num_cells for sd in subdomains])
         permeability = pp.wrap_as_dense_ad_array(
-            self.solid.fracture_permeability, size, name="fracture permeability"
+            self.solid.fracture_permeability,
+            size=size,
+            name="fracture permeability",
+            grids=subdomains,
         )
         return self.isotropic_second_order_tensor(subdomains, permeability)
 

src/porepy/examples/flow_benchmark_3d_case_2.py

@@ -121,7 +121,7 @@ def matrix_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
                 vals.append(kxx)  # This is just a placeholder, it will not be used.
 
         if len(subdomains) > 0:
-            permeability = pp.wrap_as_dense_ad_array(np.hstack(vals))
+            permeability = pp.wrap_as_dense_ad_array(np.hstack(vals), grids=subdomains)
         else:
             # This part is need it only to set up the model, not actually used.
             permeability = self.solid.permeability  # type:ignore
@@ -142,7 +142,10 @@ def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator
         # Use `fracture_permeability` as intersection permeability under the assumption
         # that they are equal. This is valid in the current benchmark case.
         permeability = pp.wrap_as_dense_ad_array(
-            self.solid.fracture_permeability, size, name="intersection_permeability"
+            self.solid.fracture_permeability,
+            size=size,
+            name="intersection_permeability",
+            grids=subdomains,
         )
         return self.isotropic_second_order_tensor(subdomains, permeability)
 

src/porepy/examples/flow_benchmark_3d_case_3.py

@@ -76,7 +76,10 @@ def intersection_permeability(self, subdomains: list[pp.Grid]) -> pp.ad.Operator
         # Use `fracture_permeability` as intersection permeability under the assumption
         # that they are equal. This is valid in the current benchmark case.
         permeability = pp.wrap_as_dense_ad_array(
-            self.solid.fracture_permeability, size, name="intersection permeability"
+            self.solid.fracture_permeability,
+            size=size,
+            name="intersection permeability",
+            grids=subdomains,
         )
         return self.isotropic_second_order_tensor(subdomains, permeability)
 

src/porepy/models/abstract_equations.py

@@ -15,12 +15,12 @@
 from __future__ import annotations
 
 from functools import cached_property
-from typing import Callable, Optional, Sequence, Union, cast
+from typing import Callable, Mapping, Optional, Sequence, Union, cast
 
 import numpy as np
 
 import porepy as pp
-from porepy.numerics.ad.equation_system import GridEntity
+from porepy.numerics.ad.equation_system import GridEntities, GridEntity
 
 
 class EquationMixin(pp.PorePyModel):
@@ -231,7 +231,9 @@ def eliminate_locally(
         dependencies: Sequence[Callable[[pp.GridLikeSequence], pp.ad.Variable]],
         func: Callable[..., tuple[np.ndarray, np.ndarray]],
         domains: Sequence[pp.Grid | pp.MortarGrid | pp.BoundaryGrid],
-        equations_per_grid_entity: Optional[dict[GridEntity, int]] = None,
+        equations_per_grid_entity: Optional[
+            Union[GridEntities, Mapping[GridEntity, int]]
+        ] = None,
     ) -> None:
         """Method to add a secondary equation eliminating a variable by some
         constitutive law depending on *other* variables.
@@ -279,12 +281,12 @@ def eliminate_locally(
             equations_per_grid_entity: ``default=None``.
 
                 Argument for when adding above equation to the equation system and
-                creating a surrogate factory. If None, the default ``{'cells':1}`` is
-                assigned.
+                creating a surrogate factory. If None, the default
+                ``GridEntities(cells=1)`` is assigned.
 
         """
         if equations_per_grid_entity is None:
-            equations_per_grid_entity = {"cells": 1}
+            equations_per_grid_entity = GridEntities(cells=1)
 
         # Separate these two because Boundary values for independent quantities are
         # stored differently.
@@ -310,11 +312,7 @@ def eliminate_locally(
         equ_name = f"elimination_of_{sec_var.name}_on_grids_{g_ids}"
         local_equ = sec_var - sec_expr(non_boundaries)
         local_equ.set_name(equ_name)
-        self.equation_system.set_equation(
-            local_equ,
-            cast(list[pp.Grid] | list[pp.MortarGrid], non_boundaries),
-            equations_per_grid_entity,
-        )
+        self.equation_system.set_equation(local_equ)
 
         self._add_local_elimination(
             equ_name, sec_var, sec_expr, func, cast(pp.GridLikeSequence, domains)

src/porepy/models/boundary_condition.py

@@ -135,7 +135,7 @@ def update_boundary_condition(
             pp.set_solution_values(name=name, values=vals, data=data, iterate_index=0)
 
     def create_boundary_operator(
-        self, name: str, domains: Sequence[pp.BoundaryGrid]
+        self, name: str, domains: Sequence[pp.BoundaryGrid], dim: int = 1
     ) -> pp.ad.TimeDependentDenseArray:
         """Creates an operator on boundary grids.
 
@@ -155,7 +155,12 @@ def create_boundary_operator(
         """
         if not all(isinstance(x, pp.BoundaryGrid) for x in domains):
             raise ValueError("Domains must consist entirely of the boundary grids.")
-        return pp.ad.TimeDependentDenseArray(name=name, domains=domains)
+        return pp.ad.TimeDependentDenseArray(
+            name=name,
+            domains=domains,
+            dof_info={pp.ad.GridEntity.cells: dim},
+            domain_type=pp.ad.DomainType.boundary_grids,
+        )
 
     def _combine_boundary_operators(
         self,
@@ -194,10 +199,16 @@ def _combine_boundary_operators(
         }
         filters = {
             "dirichlet": pp.ad.TimeDependentDenseArray(
-                name=(name + "_filter_dir"), domains=boundary_grids
+                name=(name + "_filter_dir"),
+                domains=boundary_grids,
+                dof_info={pp.ad.GridEntity.cells: dim},
+                domain_type=pp.ad.DomainType.boundary_grids,
             ),
             "neumann": pp.ad.TimeDependentDenseArray(
-                name=(name + "_filter_neu"), domains=boundary_grids
+                name=(name + "_filter_neu"),
+                domains=boundary_grids,
+                dof_info={pp.ad.GridEntity.cells: dim},
+                domain_type=pp.ad.DomainType.boundary_grids,
             ),
         }
 
@@ -206,7 +217,10 @@ def _combine_boundary_operators(
         if robin_operator is not None:
             operators["robin"] = robin_operator(boundary_grids)
             filters["robin"] = pp.ad.TimeDependentDenseArray(
-                name=(name + "_filter_rob"), domains=boundary_grids
+                name=(name + "_filter_rob"),
+                domains=boundary_grids,
+                dof_info={pp.ad.GridEntity.cells: dim},
+                domain_type=pp.ad.DomainType.boundary_grids,
             )
 
         # Adding bc_type function to local storage to evaluate it before every time step

src/porepy/models/compositional_flow.py

@@ -611,7 +611,7 @@ def set_equations(self) -> None:
         for component in self.fluid.components:
             if self.has_independent_fraction(component):
                 sd_eq = self.component_mass_balance_equation(component, subdomains)
-                self.equation_system.set_equation(sd_eq, subdomains, {"cells": 1})
+                self.equation_system.set_equation(sd_eq)
 
     def component_mass_balance_equation(
         self, component: pp.Component, subdomains: list[pp.Grid]