Line search

src/porepy/models/momentum_balance.py

@@ -790,7 +790,7 @@ def initial_condition(self) -> None:
             sd.num_cells for sd in self.mdg.subdomains(dim=self.nd - 1)
         )
         traction_vals = np.zeros((self.nd, num_frac_cells))
-        traction_vals[-1] = self.solid.convert_units(-1, "Pa")
+        traction_vals[-1] = -1  # Unitary nondimensional traction.
         self.equation_system.set_variable_values(
             traction_vals.ravel("F"),
             [self.contact_traction_variable],

src/porepy/models/solution_strategy.py

@@ -11,6 +11,7 @@
 import logging
 import time
 import warnings
+from functools import partial
 from pathlib import Path
 from typing import Any, Callable, Optional
 
@@ -724,3 +725,198 @@ def update_time_dependent_ad_arrays(self) -> None:
 
         """
         self.update_all_boundary_conditions()
+
+
+class ContactIndicators:
+    """Class for computing contact indicators used for tailored line search.
+
+    This functionality is experimental and may be subject to change.
+
+    The class is a mixin for the solution strategy classes for models with contact
+    mechanics. The class provides methods for computing the opening and sliding
+    indicators, which are used for tailored line search as defined in the class
+    :class:`~porepy.numerics.nonlinear.line_search.ConstraintLineSearch`.
+
+    By specifying the parameter `adaptive_indicator_scaling` in the model
+    parameters, the indicators can be scaled adaptively by the characteristic
+    fracture traction estimate based on the most recent iteration value.
+
+    """
+
+    basis: Callable[[list[pp.Grid], int], list[pp.ad.SparseArray]]
+    """Basis vector operator."""
+
+    contact_traction: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Contact traction operator."""
+
+    contact_mechanics_numerical_constant: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Contact mechanics numerical constant."""
+
+    displacement_jump: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Displacement jump operator."""
+
+    equation_system: pp.ad.EquationSystem
+    """Equation system manager."""
+
+    fracture_gap: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Fracture gap operator."""
+
+    friction_bound: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Friction bound operator."""
+
+    normal_component: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Normal component operator for vectors defined on fracture subdomains."""
+
+    tangential_component: Callable[[list[pp.Grid]], pp.ad.Operator]
+    """Tangential component operator for vectors defined on fracture subdomains."""
+
+    mdg: pp.MixedDimensionalGrid
+    """Mixed-dimensional grid."""
+
+    nd: int
+    """Ambient dimension of the problem."""
+
+    params: dict[str, Any]
+    """Dictionary of parameters."""
+
+    def opening_indicator(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:
+        """Function describing the state of the opening constraint.
+
+        The function is a linear combination of the two arguments of the max
+        function of the normal fracture deformation equation. Arbitrary sign
+        convention: Negative for open fractures, positive for closed ones.
+
+        The parameter `adaptive_indicator_scaling` scales the indicator by the
+        contact traction estimate.
+
+        Parameters:
+            subdomains: List of fracture subdomains.
+
+        Returns:
+            opening_indicator: Opening indicator operator.
+
+        """
+        nd_vec_to_normal = self.normal_component(subdomains)
+        # The normal component of the contact traction and the displacement jump
+        t_n: pp.ad.Operator = nd_vec_to_normal @ self.contact_traction(subdomains)
+        u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(subdomains)
+        c_num = self.contact_mechanics_numerical_constant(subdomains)
+        max_arg_1 = pp.ad.Scalar(-1.0) * t_n
+        max_arg_2 = c_num * (u_n - self.fracture_gap(subdomains))
+        ind = max_arg_1 - max_arg_2
+        if self.params.get("adaptive_indicator_scaling", False):
+            # Scale adaptively based on the contact traction estimate.
+            # Base variable values from all fracture subdomains.
+            all_subdomains = self.mdg.subdomains(dim=self.nd - 1)
+            scale_op = self.contact_traction_estimate(all_subdomains)
+            scale = self.compute_traction_norm(scale_op.value(self.equation_system))
+            ind = ind / pp.ad.Scalar(scale)
+        return ind
+
+    def sliding_indicator(
+        self,
+        subdomains: list[pp.Grid],
+    ) -> pp.ad.Operator:
+        """Function describing the state of the sliding constraint.
+
+        The function is a linear combination of the two arguments of the max function
+        of the tangential fracture deformation equation. Arbitrary sign convention:
+        Negative for sticking, positive for sliding:  ||T_t+c_t u_t||-b_p
+
+        The parameter `adaptive_indicator_scaling` scales the indicator by the contact
+        traction estimate.
+
+        Parameters:
+            subdomains: List of fracture subdomains.
+
+        Returns:
+            sliding_indicator: Sliding indicator operator.
+
+        """
+
+        # Basis vector combinations
+        num_cells = sum([sd.num_cells for sd in subdomains])
+        # Mapping from a full vector to the tangential component
+        nd_vec_to_tangential = self.tangential_component(subdomains)
+
+        tangential_basis: list[pp.ad.SparseArray] = self.basis(
+            subdomains, dim=self.nd - 1  # type: ignore[call-arg]
+        )
+
+        # Variables: The tangential component of the contact traction and the
+        # displacement jump
+        t_t: pp.ad.Operator = nd_vec_to_tangential @ self.contact_traction(subdomains)
+        u_t: pp.ad.Operator = nd_vec_to_tangential @ self.displacement_jump(subdomains)
+        # The time increment of the tangential displacement jump
+        u_t_increment: pp.ad.Operator = pp.ad.time_increment(u_t)
+        zeros_frac = pp.ad.DenseArray(np.zeros(num_cells))
+
+        f_max = pp.ad.Function(pp.ad.maximum, "max_function")
+        f_norm = pp.ad.Function(partial(pp.ad.l2_norm, self.nd - 1), "norm_function")
+        # Heaviside function. The 0 as the second argument to partial() implies
+        # f_heaviside(0)=0, a choice that is not expected to affect the result in this
+        # context.
+        f_heaviside = pp.ad.Function(partial(pp.ad.heaviside, 0), "heaviside_function")
+
+        c_num_as_scalar = self.contact_mechanics_numerical_constant(subdomains)
+
+        c_num = pp.ad.sum_operator_list(
+            [e_i * c_num_as_scalar * e_i.T for e_i in tangential_basis]
+        )
+        tangential_sum = t_t + c_num @ u_t_increment
+
+        max_arg_1 = f_norm(tangential_sum)
+        max_arg_1.set_name("norm_tangential")
+
+        max_arg_2 = f_max(self.friction_bound(subdomains), zeros_frac)
+        max_arg_2.set_name("b_p")
+
+        h_oi = f_heaviside(self.opening_indicator(subdomains))
+        ind = max_arg_1 - max_arg_2
+
+        if self.params.get("adaptive_indicator_scaling", False):
+            # Base on all fracture subdomains
+            all_subdomains = self.mdg.subdomains(dim=self.nd - 1)
+            scale_op = self.contact_traction_estimate(all_subdomains)
+            scale = self.compute_traction_norm(scale_op.value(self.equation_system))
+            ind = ind / pp.ad.Scalar(scale)
+        return ind * h_oi
+
+    def contact_traction_estimate(self, subdomains):
+        """Estimate the characteristic size of contact traction.
+
+        The characteristic fracture traction is estimated as the maximum of the
+        contact traction over the fracture subdomains.
+
+        Parameters:
+            subdomains: List of subdomains where the contact traction is defined.
+
+        Returns:
+            Characteristic fracture traction estimate.
+
+        """
+        # The normal component of the contact traction and the displacement jump
+        t: pp.ad.Operator = self.contact_traction(subdomains)
+        e_n = self.e_i(subdomains, dim=self.nd, i=self.nd - 1)
+
+        u = self.displacement_jump(subdomains) - e_n @ self.fracture_gap(subdomains)
+        c_num = self.contact_mechanics_numerical_constant(subdomains)
+        f_norm = pp.ad.Function(partial(pp.ad.l2_norm, self.nd), "norm_function")
+        return f_norm(t) + f_norm(c_num * u)
+
+    def compute_traction_norm(self, val: np.ndarray) -> float:
+        """Compute a norm of the traction estimate from the vector-valued traction.
+
+        The scalar traction is computed as the p norm of the traction vector.
+
+        Parameters:
+            val: Vector-valued traction.
+
+        Returns:
+            Scalar traction.
+
+        """
+        val = val.clip(1e-8, 1e8)
+        p = self.params.get("traction_estimate_p_mean", 5.0)
+        p_mean = np.mean(val**p, axis=0) ** (1 / p)
+        return float(p_mean)

src/porepy/numerics/ad/functions.py

@@ -254,9 +254,30 @@ def arctanh(var: FloatType) -> FloatType:
 
 
 # Step and Heaviside functions
-def heaviside(var, zerovalue: float = 0.5):
-    if isinstance(var, AdArray):
-        return np.heaviside(var.val, zerovalue)
+def heaviside(zerovalue: float, var: FloatType) -> FloatType:
+    """Heaviside function.
+
+    The Heaviside function is defined as:
+    .. math::
+        H(x) = \\begin{cases}
+            0, & x < 0, \\\\
+            zerovalue, & x = 0, \\\\
+            1, & x > 0.
+        \\end{cases}
+
+    Parameters:
+        zerovalue: Value of the Heaviside function at zero. Typically, this is
+        set to 0, 0.5 or 1.
+        var: Input array.
+
+    Returns:
+        Heaviside function (and its Jacobian if applicable) in form of a AdArray
+        or ndarray (depending on the input).
+
+    """
+    if isinstance(var, pp.ad.AdArray):
+        zero_jac = sps.csr_matrix(var.jac.shape)
+        return pp.ad.AdArray(np.heaviside(var.val, zerovalue), zero_jac)
     else:
         return np.heaviside(var, zerovalue)