MAINT: Get rid of caching md variables.

Author: jwboth

Diff for: src/porepy/models/metric.py

@@ -7,6 +7,7 @@
 
 from abc import ABC
 from functools import partial
+from typing import Callable
 
 import numpy as np
 
@@ -16,14 +17,14 @@
 class BaseMetric(ABC):
     """Base class for a metric on algebraic variant of mixed-dimensional variables."""
 
-    def variable_norm(self, solution: np.ndarray) -> float:
+    def variable_norm(self, values: np.ndarray) -> float:
         """Base method for measuring the size of physical states (increments, solution etc.)
 
         Parameters:
-            solution: algebraic respresentation of a mixed-dimensional variable
+            values: algebraic respresentation of a mixed-dimensional variable
 
         Returns:
-            float: measure of solution
+            float: measure of values
 
         """
         raise NotImplementedError("This is an abstract class.")
@@ -36,67 +37,64 @@ class EuclideanMetric(BaseMetric):
     e.g. considering errors for each variable and/or each grid separately,
     possibly using _l2_norm_cell
 
-    We normalize by the size of the solution vector as proxy for domain size.
+    We normalize by the size of the vector as proxy for domain size.
 
     """
 
-    def variable_norm(self, solution: np.ndarray) -> float:
+    def variable_norm(self, values: np.ndarray) -> float:
         """Implementation of Euclidean l2 norm of the full vector, scaled by vector size.
 
         Parameters:
-            solution: algebraic respresentation of a mixed-dimensional variable
+            values: algebraic respresentation of a mixed-dimensional variable
 
         Returns:
-            float: measure of solution
+            float: measure of values
 
         """
-        return np.linalg.norm(solution) / np.sqrt(solution.size)
+        return np.linalg.norm(values) / np.sqrt(values.size)
 
 
 class LebesgueMetric(BaseMetric):
     """Dimension-consistent Lebesgue metric (blind to physics), but separates dimensions."""
 
-    def variable_norm(self, solution: np.ndarray) -> float:
-        """Implementation of Lebesgue L2 norm of the physical solution.
+    equation_system: pp.EquationSystem
+    """EquationSystem object for the current model. Normally defined in a mixin class
+    defining the solution strategy.
+
+    """
+    volume_integral: Callable[[pp.ad.Operator, list[pp.Grid], int], pp.ad.Operator]
+    """General volume integral operator, defined in `pp.BalanceEquation`."""
+
+    def variable_norm(self, values: np.ndarray) -> float:
+        """Implementation of mixed-dimensional Lebesgue L2 norm of a physical state.
 
         Parameters:
-            solution: algebraic respresentation of a mixed-dimensional variable
+            values: algebraic respresentation of a mixed-dimensional variable
 
         Returns:
-            float: measure of solution
+            float: measure of values
 
         """
         # Initialize container for collecting separate L2 norms (squarred).
         integrals_squarred = []
 
-        # Convert algebraic solution vector to mixed dimensional variable
-        # TODO - The current implementation utilizes the datastructures in place, which is
-        # not very elegant...the creation of a isolated MixedDimensionalVariable without
-        # data transfer would be better.
-        cached_variable_values = self.equation_system.get_variable_values(
-            iterate_index=0
-        )
-        self.equation_system.set_variable_values(solution, iterate_index=0)
-
-        # Treat each variable separately - with this separate also dimensions
+        # Use the equation system to get a view onto mixed-dimensional data structures.
+        # Compute the L2 norm of each variable separately, automatically taking into
+        # account volume and specific volume
         for variable in self.equation_system.variables:
 
-            # Compute square of the L2 norm taking into account volume and specific volume
             l2_norm = pp.ad.Function(partial(pp.ad.l2_norm, variable.dim), "l2_norm")
             sd = variable.domain
+            indices = self.equation_system.dofs_of(variable)
+            ad_values = pp.ad.DenseArray(values[indices])
             integral_squarred = np.sum(
-                self.volume_integral(l2_norm(variable) ** 2, [sd], dim=1).value(
+                self.volume_integral(l2_norm(ad_values) ** 2, [sd], 1).value(
                     self.equation_system
                 )
             )
 
             # Collect the L2 norm squared.
             integrals_squarred.append(integral_squarred)
 
-        # Reset the cached variable values
-        self.equation_system.set_variable_values(
-            cached_variable_values, iterate_index=0
-        )
-
         # Squash all results by employing a consistent L2 approach.
         return np.sqrt(np.sum(integrals_squarred))