From 383a2f05ee4bc61df7a3e588e08f87a940952f39 Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Wed, 7 Aug 2024 13:16:46 +0000 Subject: [PATCH 1/8] BUG: Improve Heaviside implementation --- src/porepy/numerics/ad/functions.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/src/porepy/numerics/ad/functions.py b/src/porepy/numerics/ad/functions.py index 9f29f75f72..1df2b4a2c6 100644 --- a/src/porepy/numerics/ad/functions.py +++ b/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) From 23bf2caf35c127666c6d05c7ab48c3d6d14236db Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Wed, 7 Aug 2024 13:17:09 +0000 Subject: [PATCH 2/8] FEAT: Add line search functionality --- src/porepy/models/solution_strategy.py | 196 ++++++ src/porepy/numerics/nonlinear/line_search.py | 663 +++++++++++++++++++ 2 files changed, 859 insertions(+) create mode 100644 src/porepy/numerics/nonlinear/line_search.py diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py index 1f05154545..92e375eb73 100644 --- a/src/porepy/models/solution_strategy.py +++ b/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. + + 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. + + TODO: Consider moving the class to a more appropriate (new) module. + + """ + + 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 `relative_violation` scales the indicator by the + characteristic fracture 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 characteristic fracture traction estimate. + # Base on all fracture subdomains. + all_subdomains = self.mdg.subdomains(dim=self.nd - 1) + scale_op = self.characteristic_fracture_traction_estimate(all_subdomains) + scale = self.compute_scalar_traction(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 `relative_violation` scales the indicator by the characteristic + fracture 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 with zerovalue at argument.value=0 + 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.characteristic_fracture_traction_estimate(all_subdomains) + scale = self.compute_scalar_traction(scale_op.value(self.equation_system)) + ind = ind / pp.ad.Scalar(scale) + return ind * h_oi + + def characteristic_fracture_traction_estimate(self, subdomains): + """Estimate the characteristic fracture 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. + + """ + # The normal component of the contact traction and the displacement jump + t: pp.ad.Operator = self.contact_traction( + subdomains + ) / self.characteristic_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_scalar_traction(self, val: np.ndarray) -> float: + """Compute the scalar traction 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("characteristic_traction_p_mean", 5.0) + p_mean = np.mean(val**p, axis=0) ** (1 / p) + return float(p_mean) diff --git a/src/porepy/numerics/nonlinear/line_search.py b/src/porepy/numerics/nonlinear/line_search.py new file mode 100644 index 0000000000..dda79de158 --- /dev/null +++ b/src/porepy/numerics/nonlinear/line_search.py @@ -0,0 +1,663 @@ +"""Module for line search algorithms for nonlinear solvers. + +The algorithm is described in the manuscript at https://arxiv.org/abs/2407.01184. + +Main active classes, combined in the NonlinearSolver class: +- LineSearchNewtonSolver - extends NewtonSolver with line search and implements a basic + line search based on the residual norm. +- SplineInterpolationLineSearch - implements a line search based on spline interpolation. +- ConstraintLineSearch - implements a line search based on constraint functions for contact + mechanics. + +The functionality is invoked by specifying the solver in the model parameters, e.g.: + + ```python + model.params["nonlinear_solver"] = NonlinearSolver + ``` + +The solver can be further customized by specifying parameters in the model parameters. It also +requires implementation of the constraint functions in the model as methods called +"opening_indicator" and "sliding_indicator", see model_setup.ContactIndicators. + +""" + +import logging +from typing import Any, Callable, Optional, Sequence, cast + +import numpy as np +import scipy + +import porepy as pp + +logger = logging.getLogger(__name__) + + +class LineSearchNewtonSolver(pp.NewtonSolver): + """Class for relaxing a nonlinear iteration update using a line search. + + This class extends the iteration method to include a line search and implements + a line search based on the residual norm. + + """ + + def iteration(self, model) -> np.ndarray: + """A single nonlinear iteration. + + Add line search to the iteration method. First, call the super method to + compute a nonlinear iteration update. Then, perform a line search along + that update direction. + + Parameters: + model: The simulation model. + + Returns: + The solution update. + + """ + dx = super().iteration(model) + relaxation_vector = self.nonlinear_line_search(model, dx) + + # Update the solution + sol = relaxation_vector * dx + model._current_update = sol + return sol + + def nonlinear_line_search(self, model, dx: np.ndarray) -> np.ndarray: + """Perform a line search along the Newton step. + + Parameters: + model: The model. + dx: The nonlinear iteration update. + + Returns: + The step length vector, one for each degree of freedom. The relaxed update + is obtained by multiplying the nonlinear update step by the relaxation vector. + + """ + return self.residual_line_search(model, dx) + + def residual_line_search(self, model, dx: np.ndarray) -> np.ndarray: + """Compute the relaxation factors for the current iteration based on the residual. + + Parameters: + model: The model. + dx: The nonlinear iteration update. + + Returns: + The step length vector, one value for each degree of freedom. + + """ + if not model.params.get("Global_line_search", False): + return np.ones_like(dx) + + def objective_function(weight): + """Objective function for the trial update relaxed by + the current weight. + + Parameters: + weight: The relaxation factor. + + Returns: + The objective function value. Specified as the norm of the residual + in this implementation, may be overridden. + + """ + return self.residual_objective_function(model, dx, weight) + + tol = 1e-1 + f_0 = objective_function(0) + f_1 = objective_function(1) + if f_1 < model.params["nl_convergence_tol_res"] or (f_1 < f_0 / 1e4): + # The objective function is zero at the full nonlinear step. This + # means that the nonlinear step is a minimum of the objective + # function. We can use the update without any relaxation. + return np.ones_like(dx) + + def f_terminate(vals): + """Terminate the recursion if the objective function is increasing.""" + return vals[-1] > vals[-2] + + alpha = self.recursive_weight_from_sampling( + 0, 1, f_terminate, objective_function, f_0, f_1, 3, tol + ) + # safeguard against zero weights + return np.maximum(alpha, tol / 10) * np.ones_like(dx) + + def recursive_weight_from_sampling( + self, + a: float, + b: float, + condition_function: Callable[[Sequence], bool], + function: Callable, + f_a: Optional[float] = None, + f_b: Optional[float] = None, + num_steps: int = 5, + step_size_tolerance: float = 1e-1, + ) -> float: + """Recursive function for finding a weight satisfying a condition. + + The function is based on sampling the function in the interval [a, b] and + recursively narrowing down the interval until the interval is smaller than the + tolerance and the condition is satisfied. It returns the smallest tested value + not satisfying the condition. + + Parameters: + a, b: The interval. + condition_function: The condition function. It takes a sequence of line search + evaluations as argument and returns True if the condition is satisfied, + indicating that the recursion should be terminated. + function: The function to be tested. Returns a scalar or vector, must be compatible + with the condition function's parameter. + f_a: The value of the function at a. If not given, it is computed. + num_steps: The number of sampling points in the interval [a, b]. + step_size_tolerance: The tolerance for the step size. If the step size is smaller + than this, the recursion is terminated. + + Returns: + The smallest tested value not satisfying the condition. + + """ + if f_a is None: + f_a = function(a) + terminate_condition = False + sampling_points = np.linspace(a, b, num_steps) + step_size = (b - a) / (num_steps - 1) + f_vals = [f_a] + for c in sampling_points[1:]: + if np.isclose(c, b) and f_b is not None: + f_c = f_b + else: + f_c = function(c) + f_vals.append(f_c) + terminate_condition = condition_function(f_vals) + if not terminate_condition: + f_a = f_c + a = c + else: + # There is a local minimum in the narrowed-down interval [a, c] + if step_size > step_size_tolerance: + # Find it to better precision + return self.recursive_weight_from_sampling( + a, + c, + condition_function, + function, + f_a=f_a, + num_steps=num_steps, + step_size_tolerance=step_size_tolerance, + ) + else: + # We're happy with the precision, return the minimum + return c + + # We went through the whole interval without finding a local minimum. + # Thus, we assume that the minimum lies in [c, b]. If we have reached + # the tolerance, we return b. Otherwise, we search in [c, b]. + if step_size < step_size_tolerance: + return b + else: + return self.recursive_weight_from_sampling( + sampling_points[-2], + b, + condition_function, + function, + f_a=f_vals[-2], + num_steps=num_steps, + step_size_tolerance=step_size_tolerance, + ) + + def residual_objective_function( + self, model, dx: np.ndarray, weight: float + ) -> np.floating[Any]: + """Compute the objective function for the current iteration. + + The objective function is the norm of the residual. + + Parameters: + model: The model. + dx: The nonlinear iteration update. + weight: The relaxation factor. + + Returns: + The objective function value. + + """ + x_0 = model.equation_system.get_variable_values(iterate_index=0) + residual = model.equation_system.assemble( + state=x_0 + weight * dx, evaluate_jacobian=False + ) + return np.linalg.norm(residual) + + +class SplineInterpolationLineSearch: + """Class for computing relaxation factors based on spline interpolation. + + This class could be seen as a tool for the technical step of performing a + line search. It also specifies that this choice should be used for the + constraint weights. + + """ + + def compute_constraint_violation_weights( + self, + model, + solution_update: np.ndarray, + constraint_function: pp.ad.Operator, + crossing_inds: np.ndarray, + f_0: np.ndarray, + max_weight: float = 1.0, + interval_target_size=1e-3, + ) -> float: + """Specify that the constraint weights are computed based on spline interpolation. + + Parameters: + model: The model. + solution_update: The nonlinear iteration update. + constraint_function: The constraint function. Returns a scalar or vector and + changes sign at the zero crossing. + crossing_inds: The indices where the constraint function has changed sign. + f_0: The value of the constraint function at the initial point. + max_weight: The maximum value for the constraint weights. + interval_target_size: The target size of the interval for the root finding. The + recursion is terminated when the interval is smaller than this. + + Returns: + The step length weight. + + """ + if not np.any(crossing_inds): + return 1.0 + # If the indicator has changed, we need to compute the relaxation + # factors. We do this by recursively narrowing down the interval until + # the interval is smaller than the tolerance using a spline + # interpolation. + a, b = 0.0, max_weight + x_0 = model.equation_system.get_variable_values(iterate_index=0) + f_0 = f_0[crossing_inds] + # Compute the constraint function at the new weights. Specify that return value + # is an array. + f_1_vals = cast( + np.ndarray, + constraint_function.value(model.equation_system, x_0 + solution_update * b), + ) + f_1 = f_1_vals[crossing_inds] + + def f(x): + return constraint_function.value( + model.equation_system, x_0 + solution_update * x + )[crossing_inds] + + alpha, a, b = self.recursive_spline_interpolation( + a, + b, + f, + f_0, + f_1, + interval_target_size=interval_target_size, + method="roots", + ) + return alpha + + def recursive_spline_interpolation( + self, + a: float, + b: float, + function: Callable, + f_a: Optional[float | np.ndarray] = None, + f_b: Optional[float | np.ndarray] = None, + num_pts: int = 5, + interval_target_size: float = 1e-1, + method="minimize_scalar", + ) -> tuple[float, float, float]: + """Recursive function for finding a weight satisfying a condition. + + Returns both the optimum/root (see method) and the minimal interval in + which it is assumed to lie. + + Parameters: + a, b: The interval. + function: The function to be tested. Returns a scalar or vector defined on + the interval [a, b]. + f_a: The value of the function at a. If not given, it is computed. + f_b: The value of the function at b. If not given, it is computed. + num_pts: The number of sampling points in the interval [a, b]. + step_size_tolerance: The tolerance for the step size. If the step size is smaller + than this, the recursion is terminated. + method: The method for finding the minimum of the spline. Either "minimize_scalar" + or "roots". + + Returns: + Tuple containing: + The minimum of the function in the interval [a, b]. + The lower bound of the interval in which the minimum is assumed to lie. + The upper bound of the interval in which the minimum is assumed to lie. + """ + counter = 0 + while b - a > interval_target_size or counter < 1: + alpha, x, y = self.optimum_from_spline( + function, + a, + b, + f_a, + f_b, + num_pts=num_pts, + method=method, + ) + x = np.linspace(a, b, num_pts) + # Find the indices on either side of alpha + ind = np.searchsorted(x, alpha) + if ind == 0: + b = x[1] + f_b = y[1] + elif ind == num_pts: + a = x[ind - 1] + f_a = y[ind - 1] + else: + a = x[ind - 1] + b = x[ind] + f_a = y[ind - 1] + f_b = y[ind] + counter += 1 + + return alpha, a, b + + def optimum_from_spline( + self, + f: Callable, + a: float, + b: float, + f_a=None, + f_b=None, + num_pts: int = 5, + method="minimize_scalar", + ) -> tuple[float, np.ndarray, np.ndarray]: + """Compute the minimum/root of the spline interpolation of the function. + + Parameters: + f: The function to be interpolated. + a, b: The interval. + f_a: The value of the function at a. If not given, it is computed. + f_b: The value of the function at b. If not given, it is computed. + num_pts: The number of sampling points in the interval [a, b]. + method: The method for finding the minimum of the spline. Either "minimize_scalar" + or "roots". + + Returns: + Tuple containing: + The minimum of the function in the interval [a, b]. + The x values of the spline interpolation. + The y values of the spline interpolation. + + """ + + x = np.linspace(a, b, num_pts) + y_list = [] + + for pt in x: + if f_a is not None and np.isclose(pt, a): + f_pt = f_a + elif f_b is not None and np.isclose(pt, b): + f_pt = f_b + else: + f_pt = f(pt) + if np.any(np.isnan(f_pt)): + # If we get overflow, truncate the x vector + x = x[: np.where(x == pt)[0][0]] + break + # Collect function values, scalar or vector + y_list.append(f_pt) + if isinstance(y_list[0], np.ndarray): + y = np.vstack(y_list) + else: + y = np.array(y_list) + + def compute_and_postprocess_single(poly, a: float, b: float) -> float: + if method == "minimize_scalar": + minimum = scipy.optimize.minimize_scalar( + lambda s: poly(s), bounds=[a, b], method="bounded" + ) + min_x = minimum.x + elif method == "roots": + min_x = poly.roots() + if min_x.size == 0: + return b + else: + # Find smallest root inside [a, b] + min_x = min_x[(min_x >= a) & (min_x <= b)] + if min_x.size == 0: + return b + else: + return np.min(min_x) + + # Find minima of the spline + + if isinstance(y_list[0], np.ndarray): + all_minima = [] + for i in range(y.shape[1]): + poly = scipy.interpolate.PchipInterpolator(x, y[:, i]) + this_min = compute_and_postprocess_single(poly, a, b) + all_minima.append(this_min) + alpha = np.min(all_minima) + else: + poly = scipy.interpolate.PchipInterpolator(x, y) + alpha = compute_and_postprocess_single(poly, a, b) + + return alpha, x, y + + +class ConstraintLineSearch: + """Class for computing relaxation weights based on constraint functions + for contact mechanics. + + The contract with the Model class is that the constraint functions are + defined in the model as Operator returning methods called "opening_indicator" + and "sliding_indicator". + + """ + + compute_constraint_violation_weights: Callable[..., float] + """Method for computing the constraint weights. + + This method specifies the algorithm for computing the constraint weights by a line + search method. Current option is based on spline interpolation. + + """ + + residual_line_search: Callable[[Any, np.ndarray], np.ndarray] + """Method for computing the relaxation factors for the current iteration + based on the residual.""" + + @property + def use_fracture_minimum(self): + return True + + @property + def min_line_search_weight(self): + """Minimum weight for the relaxation weights.""" + return 1e-12 + + def nonlinear_line_search(self, model, dx: np.ndarray) -> np.ndarray: + """Perform a line search along the Newton step. + + First, call super method using the global residual as the objective function. + Then, compute the constraint weights. + + Parameters: + model: The model. + dx: The Newton step. + + Returns: + The step length vector, one for each degree of freedom. + + """ + residual_weight = self.residual_line_search(model, dx) + if model.params.get("Local_line_search", False): + return self.constraint_line_search(model, dx, residual_weight.min()) + else: + return residual_weight + + def constraint_line_search( + self, model, dx: np.ndarray, max_weight: float = 1 + ) -> np.ndarray: + """Perform line search along the Newton step for the constraints. + + This method defines the constraint weights for the contact mechanics and + how they are combined to a global weight. For more advanced combinations, + this method can be overridden or the stored weights can be accessed elsewhere. + + Parameters: + model: The model. + dx: The Newton step. + max_weight: The maximum weight for the constraint weights. + + Returns: + The step length vector, one for each degree of freedom. + + """ + + subdomains = model.mdg.subdomains(dim=model.nd - 1) + + sd_weights = {} + global_weight = max_weight + for sd in subdomains: + sd_list = [sd] + # Compute the relaxation factors for the normal and tangential + # component of the contact mechanics. + normal_weights = self.constraint_weights( + model, + dx, + model.opening_indicator(sd_list), + max_weight=max_weight, + ) + tangential_weights = self.constraint_weights( + model, + dx, + model.sliding_indicator(sd_list), + max_weight=np.minimum(max_weight, normal_weights).min(), + ) + # For each cell, take minimum of tangential and normal weights + combined_weights = np.vstack((tangential_weights, normal_weights)) + min_weights = np.min(combined_weights, axis=0) + # Store the weights for the subdomain. This facilitates more advanced + # combinations of the constraint weights, e.g. using local weights for + # each cell. + sd_weights[sd] = min_weights + model.mdg.subdomain_data(sd).update({"constraint_weights": min_weights}) + # Combine the weights for all subdomains to a global minimum weight. + global_weight = np.minimum(global_weight, min_weights.min()) + # Return minimum of all weights. + weight = np.ones_like(dx) * global_weight + return weight + + def constraint_weights( + self, + model, + solution_update: np.ndarray, + constraint_function: pp.ad.Operator, + max_weight: float = 1, + ) -> np.ndarray: + """Compute weights for a given constraint. + + This method specifies the algorithm for computing the constraint weights: + - Identify the indices where the constraint function has changed sign. + - Compute the relaxation factors for these indices, allowing transition beyond + zero by a tolerance given by the constraint_violation_tolerance parameter. + - Reassess the constraint function at the new weights and repeat the process if + too many indices are still transitioning. + + Parameters: + model: The model. + solution_update: The Newton step. + constraint_function: The constraint function. + max_weight: The maximum weight for the constraint weights. + + Returns: + The step length vector, one for each degree of freedom of the constraint. + + """ + # If the sign of the function defining the regions has not changed, we + # use unitary relaxation factors + x_0 = model.equation_system.get_variable_values(iterate_index=0) + violation_tol = model.params.get("constraint_violation_tolerance", 3e-1) + relative_cell_tol = model.params.get( + "relative_constraint_transition_tolerance", 2e-1 + ) + # Compute the constraint function at the maximum weights. Specify that return value + # is an array to avoid type errors. + f_1 = cast( + np.ndarray, + constraint_function.value( + model.equation_system, x_0 + max_weight * solution_update + ), + ) + weight = max_weight + weights = max_weight * np.ones(f_1.shape) + f_0 = constraint_function.value(model.equation_system, x_0) + active_inds = np.ones(f_1.shape, dtype=bool) + for i in range(10): + # Only consider dofs where the constraint violation has changed sign + violation = violation_tol * np.sign(f_1) + f = constraint_function - pp.wrap_as_dense_ad_array(violation) + # Absolute tolerance should be safe, as constraints are assumed to be + # scaled to approximately 1. + roundoff = 1e-8 + inds = np.logical_and(np.abs(f_1) > violation_tol, f_0 * f_1 < -roundoff) + if i > 0: # Ensure at least one iteration. + if sum(active_inds) < max(1, relative_cell_tol * active_inds.size): + # Only a few indices are active, and the set of active indices + # does not contain any new indices. We can stop the iteration. + break + + else: + logger.info( + f"Relaxation factor {weight} is too large for {sum(active_inds)}" + + " indices. Reducing constraint violation tolerance." + ) + + f_0_v = f_0 - violation + # Compute the relaxation factors for the indices where the constraint + # violation has changed sign. The weight corresponds to the point at + # which the constraint function crosses zero. + crossing_weight = self.compute_constraint_violation_weights( + model, + solution_update, + f, + inds, + f_0_v, + max_weight=max_weight, + interval_target_size=1e-3, + ) + weight = np.clip( + crossing_weight, a_max=max_weight, a_min=self.min_line_search_weight + ) + weights[inds] = weight + + if not self.use_fracture_minimum: + break # Experimental. + # Check how many indices are active for current weight + f_1 = cast( + np.ndarray, + constraint_function.value( + model.equation_system, x_0 + weight * solution_update + ), + ) + active_inds = np.logical_and( + np.abs(f_1) > violation_tol, f_0 * f_1 < -roundoff + ) + if i == 9: + logger.info( + "Maximum number of iterations reached. " + + "Returning current weights." + ) + max_weight = weight + violation_tol = violation_tol / 2 + + return weights + + +# class NonlinearSolver( +# ConstraintLineSearch, +# SplineInterpolationLineSearch, +# LineSearchNewtonSolver, +# ): +# pass From 34ff3f1529540a95b249d9d1b1303d033b92dc08 Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Wed, 7 Aug 2024 16:08:36 +0000 Subject: [PATCH 3/8] BUG: Wrong method name. --- src/porepy/models/solution_strategy.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py index 92e375eb73..b3d6381dd4 100644 --- a/src/porepy/models/solution_strategy.py +++ b/src/porepy/models/solution_strategy.py @@ -896,7 +896,7 @@ def characteristic_fracture_traction_estimate(self, subdomains): # The normal component of the contact traction and the displacement jump t: pp.ad.Operator = self.contact_traction( subdomains - ) / self.characteristic_traction(subdomains) + ) / self.characteristic_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) @@ -917,6 +917,6 @@ def compute_scalar_traction(self, val: np.ndarray) -> float: """ val = val.clip(1e-8, 1e8) - p = self.params.get("characteristic_traction_p_mean", 5.0) + p = self.params.get("traction_estimate_p_mean", 5.0) p_mean = np.mean(val**p, axis=0) ** (1 / p) return float(p_mean) From 36d852b31f67d645fce82df709340598704d17bd Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Wed, 7 Aug 2024 16:09:03 +0000 Subject: [PATCH 4/8] TUT: Add tutorial on solution strategies --- tutorials/solution_strategies.ipynb | 224 ++++++++++++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100644 tutorials/solution_strategies.ipynb diff --git a/tutorials/solution_strategies.ipynb b/tutorials/solution_strategies.ipynb new file mode 100644 index 0000000000..7072a60038 --- /dev/null +++ b/tutorials/solution_strategies.ipynb @@ -0,0 +1,224 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solution strategies\n", + "This tutorials demonstrates some of the more advanced solution strategies available in `PorePy`. \n", + "The model problems are presented succinctly.\n", + "\n", + "## Line search for multiphysics problems with fracture deformation\n", + "The combination of highly nonlinear multiphysics problems and the non-smooth formulation of contact mechanics may lead to severe convergence issues.\n", + "To improve convergence, we provide a line search algorithm tailored to the irregularities of the contact equations, see https://arxiv.org/abs/2407.01184.\n", + "We demonstrate how to invoke the functionality to achieve convergence in a simple case which defies convergence for a regular Newton algorithm.\n", + "\n", + "To define a suitably challenging problem containing open, sticking and slipping fracture cells, we prescribe heterogeneous Dirichlet values for displacement on a unit cube domain containing a single fracture. \n", + "We also material parameters for granite and water, with bespoke (nonzero) values for fracture deformation parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "import porepy as pp\n", + "import numpy as np\n", + "\n", + "from porepy.applications.md_grids.model_geometries import CubeDomainOrthogonalFractures\n", + "from porepy.applications.boundary_conditions.model_boundary_conditions import BoundaryConditionsMechanicsDirNorthSouth\n", + "from porepy.numerics.nonlinear import line_search\n", + "\n", + "class DisplacementBoundaryConditions(BoundaryConditionsMechanicsDirNorthSouth):\n", + " time_manager: pp.TimeManager\n", + "\n", + " def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:\n", + " \"\"\"Boundary values for mechanics.\n", + "\n", + " Parameters:\n", + " subdomains: List of subdomains on which to define boundary conditions.\n", + "\n", + " Returns:\n", + " Array of boundary values.\n", + "\n", + " \"\"\"\n", + "\n", + " # Default is zero.\n", + " vals = np.zeros((self.nd, boundary_grid.num_cells))\n", + " if boundary_grid.dim < self.nd - 1:\n", + " return vals.ravel(\"F\")\n", + " boundary_faces = self.domain_boundary_sides(boundary_grid)\n", + "\n", + " faces = boundary_faces.north\n", + " # Start by offsetting by the fracture gap.\n", + " vals[1, faces] = self.solid.fracture_gap()\n", + " if self.time_manager.time < 1e-5:\n", + " return vals.ravel(\"F\")\n", + "\n", + " u_char = self.characteristic_displacement([boundary_grid]).value(\n", + " self.equation_system\n", + " )\n", + " # Normal component of boundary displacement.\n", + " # Produce different contact regimes along the fracture.\n", + " coo = (\n", + " boundary_grid.cell_centers[0, faces]\n", + " + boundary_grid.cell_centers[2, faces]\n", + " ) / np.max(self.domain.side_lengths())\n", + " linear_increase = u_char * (coo - 0.5)\n", + " offset = -0.5 * u_char\n", + " vals[1, faces] += offset - 2.0 * linear_increase\n", + " # Tangential component of boundary displacement.\n", + " w = self.params.get(\"tangential_u_weight\", 2)\n", + " vals[0, faces] = w * u_char * 2.0\n", + " vals[1, faces] = w * u_char * 1.0\n", + " return vals.ravel(\"F\")\n", + "\n", + "class NonzeroInitialCondition:\n", + " \"\"\"A decent initial condition for the contact problem greatly improves convergence.\"\"\"\n", + " def initial_condition(self) -> None:\n", + " \"\"\"Set the initial condition for the problem.\"\"\"\n", + " super().initial_condition()\n", + " for var in self.equation_system.variables:\n", + " if hasattr(self, \"initial_\" + var.name):\n", + " values = getattr(self, \"initial_\" + var.name)([var.domain])\n", + " self.equation_system.set_variable_values(\n", + " values, [var], iterate_index=0, time_step_index=0\n", + " )\n", + "\n", + " def initial_t(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:\n", + " \"\"\"Initial contact traction [Pa].\n", + "\n", + " Parameters:\n", + " subdomains: List of subdomains.\n", + "\n", + " Returns:\n", + " Operator for initial contact traction.\n", + "\n", + " \"\"\"\n", + " sd = subdomains[0]\n", + " traction_vals = np.zeros((self.nd, sd.num_cells))\n", + " traction_vals[-1] = -self.characteristic_contact_traction(subdomains).value(\n", + " self.equation_system\n", + " ) / self.params.get(\"characteristic_displacement_scaling\", 1.0)\n", + " return traction_vals.ravel(\"F\")\n", + "\n", + "\n", + "\n", + "\n", + "class TailoredPoromechanics(\n", + " DisplacementBoundaryConditions,\n", + " CubeDomainOrthogonalFractures,\n", + " pp.constitutive_laws.CubicLawPermeability,\n", + " pp.models.solution_strategy.ContactIndicators,\n", + " NonzeroInitialCondition,\n", + " pp.poromechanics.Poromechanics,\n", + "):\n", + " \"\"\"Combine mixins with the poromechanics class.\"\"\"\n", + "\n", + "class NonlinearSolver(\n", + " line_search.ConstraintLineSearch,\n", + " line_search.SplineInterpolationLineSearch,\n", + " line_search.LineSearchNewtonSolver,\n", + "):\n", + " \"\"\"Collect all the line search methods in one class.\"\"\"\n", + "\n", + "granite_values: dict[str, float] = pp.solid_values.granite\n", + "granite_values.update(\n", + " {\n", + " \"permeability\": 2e-15,\n", + " \"normal_permeability\": 1e-7,\n", + " \"residual_aperture\": 2e-5,\n", + " \"dilation_angle\": 0.2,\n", + " \"characteristic_displacement\": 1e-2,\n", + " }\n", + ")\n", + "params = {\n", + " \"times_to_export\": [], # No output\n", + " \"material_constants\": {\"solid\": pp.SolidConstants(granite_values), \"fluid\": pp.FluidConstants(pp.fluid_values.water)},\n", + " \"fracture_indices\": [1], # Fracture with constant y-coordinate\n", + " \"meshing_arguments\": {\"cell_size\": 1 / 4},\n", + " # Specify iterative solver parameters\n", + " \"max_iterations\": 30,\n", + " \"nl_convergence_tol\": 1e-7, # Limited by condition number\n", + " \"nl_convergence_tol_res\": 1e-7,\n", + " \"linear_solver\": \"scipy_sparse\",\n", + "}\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To invoke the line search, we need to specify a few parameters, including the \"nonlinear_solver\". \n", + "We do this in a copy of the parameter dictionary, allowing us to run two simulations." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged in 18 iterations\n" + ] + } + ], + "source": [ + "params.update(\n", + " {\n", + " \"nonlinear_solver\": NonlinearSolver,\n", + " \"Local_line_search\": 1, # Set to 0 to use turn off the tailored line search\n", + " \"Global_line_search\": 0, # Set to 1 to use turn on a residual-based line search\n", + " \"adaptive_indicator_scaling\": 1,\n", + " \"time_manager\": pp.TimeManager(\n", + " [0.0, 1e-3], # Time interval\n", + " 1e-3, # Time step\n", + " True,\n", + " ),\n", + " }\n", + ")\n", + "\n", + "model = TailoredPoromechanics(params)\n", + "pp.run_time_dependent_model(model, params)\n", + "print(\"Converged in {} iterations\".format(model.nonlinear_solver_statistics.num_iteration))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some known issues to consider:\n", + "- The nonlinearity caused by the fracture permeability has proven sensitive to discretization method.\n", + "TPFA seems considerably more stable than MPFA.\n", + "For details, see the above mentioned paper.\n", + "- If regular Newton converges, it may be faster than the line search method.\n", + "An unexplored option is to use the line search as a fallback if Newton fails to converge." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 105cd4fe584d1243afa55d1332bc06ca9d76951c Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Fri, 9 Aug 2024 09:46:06 +0000 Subject: [PATCH 5/8] TUT: Simplify line search tutorial --- tutorials/solution_strategies.ipynb | 166 ++++++---------------------- 1 file changed, 34 insertions(+), 132 deletions(-) diff --git a/tutorials/solution_strategies.ipynb b/tutorials/solution_strategies.ipynb index 7072a60038..59a180128f 100644 --- a/tutorials/solution_strategies.ipynb +++ b/tutorials/solution_strategies.ipynb @@ -6,143 +6,47 @@ "source": [ "# Solution strategies\n", "This tutorials demonstrates some of the more advanced solution strategies available in `PorePy`. \n", + "The tutorial will be supplemented with more strategies as these become available.\n", "The model problems are presented succinctly.\n", "\n", - "## Line search for multiphysics problems with fracture deformation\n", + "## Line search for multiphysics problems with fracture deformation (experimental)\n", "The combination of highly nonlinear multiphysics problems and the non-smooth formulation of contact mechanics may lead to severe convergence issues.\n", "To improve convergence, we provide a line search algorithm tailored to the irregularities of the contact equations, see https://arxiv.org/abs/2407.01184.\n", - "We demonstrate how to invoke the functionality to achieve convergence in a simple case which defies convergence for a regular Newton algorithm.\n", "\n", - "To define a suitably challenging problem containing open, sticking and slipping fracture cells, we prescribe heterogeneous Dirichlet values for displacement on a unit cube domain containing a single fracture. \n", - "We also material parameters for granite and water, with bespoke (nonzero) values for fracture deformation parameters." + "\n", + "We define a dummy problem with one fracture and trivial problem data (boundary conditions etc.). \n", + "This allows us to emphasize the steps needed to invoke the line search algorithm.\n", + "For more challenging problems which would not converge with regular Newton, please confer the examples of the above manuscript." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import porepy as pp\n", - "import numpy as np\n", "\n", "from porepy.applications.md_grids.model_geometries import CubeDomainOrthogonalFractures\n", - "from porepy.applications.boundary_conditions.model_boundary_conditions import BoundaryConditionsMechanicsDirNorthSouth\n", "from porepy.numerics.nonlinear import line_search\n", "\n", - "class DisplacementBoundaryConditions(BoundaryConditionsMechanicsDirNorthSouth):\n", - " time_manager: pp.TimeManager\n", - "\n", - " def bc_values_displacement(self, boundary_grid: pp.BoundaryGrid) -> np.ndarray:\n", - " \"\"\"Boundary values for mechanics.\n", - "\n", - " Parameters:\n", - " subdomains: List of subdomains on which to define boundary conditions.\n", - "\n", - " Returns:\n", - " Array of boundary values.\n", - "\n", - " \"\"\"\n", - "\n", - " # Default is zero.\n", - " vals = np.zeros((self.nd, boundary_grid.num_cells))\n", - " if boundary_grid.dim < self.nd - 1:\n", - " return vals.ravel(\"F\")\n", - " boundary_faces = self.domain_boundary_sides(boundary_grid)\n", - "\n", - " faces = boundary_faces.north\n", - " # Start by offsetting by the fracture gap.\n", - " vals[1, faces] = self.solid.fracture_gap()\n", - " if self.time_manager.time < 1e-5:\n", - " return vals.ravel(\"F\")\n", - "\n", - " u_char = self.characteristic_displacement([boundary_grid]).value(\n", - " self.equation_system\n", - " )\n", - " # Normal component of boundary displacement.\n", - " # Produce different contact regimes along the fracture.\n", - " coo = (\n", - " boundary_grid.cell_centers[0, faces]\n", - " + boundary_grid.cell_centers[2, faces]\n", - " ) / np.max(self.domain.side_lengths())\n", - " linear_increase = u_char * (coo - 0.5)\n", - " offset = -0.5 * u_char\n", - " vals[1, faces] += offset - 2.0 * linear_increase\n", - " # Tangential component of boundary displacement.\n", - " w = self.params.get(\"tangential_u_weight\", 2)\n", - " vals[0, faces] = w * u_char * 2.0\n", - " vals[1, faces] = w * u_char * 1.0\n", - " return vals.ravel(\"F\")\n", - "\n", - "class NonzeroInitialCondition:\n", - " \"\"\"A decent initial condition for the contact problem greatly improves convergence.\"\"\"\n", - " def initial_condition(self) -> None:\n", - " \"\"\"Set the initial condition for the problem.\"\"\"\n", - " super().initial_condition()\n", - " for var in self.equation_system.variables:\n", - " if hasattr(self, \"initial_\" + var.name):\n", - " values = getattr(self, \"initial_\" + var.name)([var.domain])\n", - " self.equation_system.set_variable_values(\n", - " values, [var], iterate_index=0, time_step_index=0\n", - " )\n", - "\n", - " def initial_t(self, subdomains: list[pp.Grid]) -> pp.ad.Operator:\n", - " \"\"\"Initial contact traction [Pa].\n", - "\n", - " Parameters:\n", - " subdomains: List of subdomains.\n", - "\n", - " Returns:\n", - " Operator for initial contact traction.\n", - "\n", - " \"\"\"\n", - " sd = subdomains[0]\n", - " traction_vals = np.zeros((self.nd, sd.num_cells))\n", - " traction_vals[-1] = -self.characteristic_contact_traction(subdomains).value(\n", - " self.equation_system\n", - " ) / self.params.get(\"characteristic_displacement_scaling\", 1.0)\n", - " return traction_vals.ravel(\"F\")\n", - "\n", - "\n", - "\n", "\n", "class TailoredPoromechanics(\n", - " DisplacementBoundaryConditions,\n", + " # Problem definition\n", " CubeDomainOrthogonalFractures,\n", + " # Add nonlinearity to the fracture flow\n", " pp.constitutive_laws.CubicLawPermeability,\n", + " # Needed for the tailored line search algorithm\n", " pp.models.solution_strategy.ContactIndicators,\n", - " NonzeroInitialCondition,\n", + " # Base class\n", " pp.poromechanics.Poromechanics,\n", "):\n", " \"\"\"Combine mixins with the poromechanics class.\"\"\"\n", "\n", - "class NonlinearSolver(\n", - " line_search.ConstraintLineSearch,\n", - " line_search.SplineInterpolationLineSearch,\n", - " line_search.LineSearchNewtonSolver,\n", - "):\n", - " \"\"\"Collect all the line search methods in one class.\"\"\"\n", "\n", - "granite_values: dict[str, float] = pp.solid_values.granite\n", - "granite_values.update(\n", - " {\n", - " \"permeability\": 2e-15,\n", - " \"normal_permeability\": 1e-7,\n", - " \"residual_aperture\": 2e-5,\n", - " \"dilation_angle\": 0.2,\n", - " \"characteristic_displacement\": 1e-2,\n", - " }\n", - ")\n", "params = {\n", - " \"times_to_export\": [], # No output\n", - " \"material_constants\": {\"solid\": pp.SolidConstants(granite_values), \"fluid\": pp.FluidConstants(pp.fluid_values.water)},\n", " \"fracture_indices\": [1], # Fracture with constant y-coordinate\n", " \"meshing_arguments\": {\"cell_size\": 1 / 4},\n", - " # Specify iterative solver parameters\n", - " \"max_iterations\": 30,\n", - " \"nl_convergence_tol\": 1e-7, # Limited by condition number\n", - " \"nl_convergence_tol_res\": 1e-7,\n", - " \"linear_solver\": \"scipy_sparse\",\n", "}\n" ] }, @@ -150,41 +54,34 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To invoke the line search, we need to specify a few parameters, including the \"nonlinear_solver\". \n", - "We do this in a copy of the parameter dictionary, allowing us to run two simulations." + "To invoke the line search, we first define the tailored Newton solver class. We then specify a few parameters, including the \"nonlinear_solver\". " ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 10, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Converged in 18 iterations\n" - ] - } - ], + "outputs": [], "source": [ + "class ConstraintLineSearchNonlinearSolver(\n", + " line_search.ConstraintLineSearch, # The tailoring to contact constraints.\n", + " line_search.SplineInterpolationLineSearch, # Technical implementation of the actual search along given update direction\n", + " line_search.LineSearchNewtonSolver, # General line search.\n", + "):\n", + " \"\"\"Collect all the line search methods in one class.\"\"\"\n", + "\n", "params.update(\n", " {\n", - " \"nonlinear_solver\": NonlinearSolver,\n", - " \"Local_line_search\": 1, # Set to 0 to use turn off the tailored line search\n", - " \"Global_line_search\": 0, # Set to 1 to use turn on a residual-based line search\n", - " \"adaptive_indicator_scaling\": 1,\n", - " \"time_manager\": pp.TimeManager(\n", - " [0.0, 1e-3], # Time interval\n", - " 1e-3, # Time step\n", - " True,\n", - " ),\n", - " }\n", + " \"nl_convergence_tol_res\": 1e-10, # Since the line search reduces the update, a residual based tolerance is needed\n", + " \"nonlinear_solver\": ConstraintLineSearchNonlinearSolver,\n", + " \"Local_line_search\": 1, # Set to 0 to use turn off the tailored line search\n", + " \"Global_line_search\": 0, # Set to 1 to use turn on a residual-based line search\n", + " \"adaptive_indicator_scaling\": 1, # Scale the indicator adaptively to increase robustness\n", + " }\n", ")\n", - "\n", + "# Run the simulation.\n", "model = TailoredPoromechanics(params)\n", - "pp.run_time_dependent_model(model, params)\n", - "print(\"Converged in {} iterations\".format(model.nonlinear_solver_statistics.num_iteration))" + "pp.run_time_dependent_model(model, params)" ] }, { @@ -198,6 +95,11 @@ "- If regular Newton converges, it may be faster than the line search method.\n", "An unexplored option is to use the line search as a fallback if Newton fails to converge." ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] } ], "metadata": { From 72929ea7cf1d0a7d824ecd3a2ad1f8fbb4739729 Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Fri, 9 Aug 2024 09:48:13 +0000 Subject: [PATCH 6/8] BUG: Remnants from code migration --- src/porepy/models/momentum_balance.py | 2 +- src/porepy/models/solution_strategy.py | 41 +++++++++++++------------- 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/porepy/models/momentum_balance.py b/src/porepy/models/momentum_balance.py index 342762f257..d8ff54e07a 100644 --- a/src/porepy/models/momentum_balance.py +++ b/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], diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py index b3d6381dd4..b426c05d59 100644 --- a/src/porepy/models/solution_strategy.py +++ b/src/porepy/models/solution_strategy.py @@ -730,6 +730,8 @@ def update_time_dependent_ad_arrays(self) -> None: 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 @@ -739,8 +741,6 @@ class ContactIndicators: parameters, the indicators can be scaled adaptively by the characteristic fracture traction estimate based on the most recent iteration value. - TODO: Consider moving the class to a more appropriate (new) module. - """ basis: Callable[[list[pp.Grid], int], list[pp.ad.SparseArray]] @@ -786,8 +786,8 @@ def opening_indicator(self, subdomains: list[pp.Grid]) -> pp.ad.Operator: function of the normal fracture deformation equation. Arbitrary sign convention: Negative for open fractures, positive for closed ones. - The parameter `relative_violation` scales the indicator by the - characteristic fracture traction estimate. + The parameter `adaptive_indicator_scaling` scales the indicator by the + contact traction estimate. Parameters: subdomains: List of fracture subdomains. @@ -805,11 +805,11 @@ def opening_indicator(self, subdomains: list[pp.Grid]) -> pp.ad.Operator: 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 characteristic fracture traction estimate. - # Base on all fracture subdomains. + # 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.characteristic_fracture_traction_estimate(all_subdomains) - scale = self.compute_scalar_traction(scale_op.value(self.equation_system)) + 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 @@ -823,8 +823,8 @@ def sliding_indicator( 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 `relative_violation` scales the indicator by the characteristic - fracture traction estimate. + The parameter `adaptive_indicator_scaling` scales the indicator by the contact + traction estimate. Parameters: subdomains: List of fracture subdomains. @@ -853,7 +853,8 @@ def sliding_indicator( 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 with zerovalue at argument.value=0 + # 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) @@ -875,13 +876,13 @@ def sliding_indicator( 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.characteristic_fracture_traction_estimate(all_subdomains) - scale = self.compute_scalar_traction(scale_op.value(self.equation_system)) + 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 characteristic_fracture_traction_estimate(self, subdomains): - """Estimate the characteristic fracture traction. + 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. @@ -890,13 +891,11 @@ def characteristic_fracture_traction_estimate(self, subdomains): subdomains: List of subdomains where the contact traction is defined. Returns: - Characteristic fracture traction. + Characteristic fracture traction estimate. """ # The normal component of the contact traction and the displacement jump - t: pp.ad.Operator = self.contact_traction( - subdomains - ) / self.characteristic_contact_traction(subdomains) + 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) @@ -904,8 +903,8 @@ def characteristic_fracture_traction_estimate(self, 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_scalar_traction(self, val: np.ndarray) -> float: - """Compute the scalar traction from the vector-valued traction. + 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. From ac58cfec5caa43974ad27348e94a03e7b018ea98 Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Fri, 9 Aug 2024 09:48:57 +0000 Subject: [PATCH 7/8] DOC: Formatting and minor improvements --- src/porepy/numerics/nonlinear/line_search.py | 198 ++++++++++--------- 1 file changed, 100 insertions(+), 98 deletions(-) diff --git a/src/porepy/numerics/nonlinear/line_search.py b/src/porepy/numerics/nonlinear/line_search.py index dda79de158..00ce049a80 100644 --- a/src/porepy/numerics/nonlinear/line_search.py +++ b/src/porepy/numerics/nonlinear/line_search.py @@ -1,23 +1,30 @@ -"""Module for line search algorithms for nonlinear solvers. +"""Module for line search algorithms for nonlinear solvers. Experimental! The algorithm is described in the manuscript at https://arxiv.org/abs/2407.01184. Main active classes, combined in the NonlinearSolver class: - LineSearchNewtonSolver - extends NewtonSolver with line search and implements a basic line search based on the residual norm. -- SplineInterpolationLineSearch - implements a line search based on spline interpolation. -- ConstraintLineSearch - implements a line search based on constraint functions for contact - mechanics. +- SplineInterpolationLineSearch - implements a line search using spline interpolation. +- ConstraintLineSearch - implements a line search based on constraint functions for + contact mechanics. The functionality is invoked by specifying the solver in the model parameters, e.g.: ```python - model.params["nonlinear_solver"] = NonlinearSolver + class ConstraintLineSearchNonlinearSolver( + ConstraintLineSearch, + SplineInterpolationLineSearch, + LineSearchNewtonSolver, + ): + pass + model.params["nonlinear_solver"] = ConstraintLineSearchNonlinearSolver ``` -The solver can be further customized by specifying parameters in the model parameters. It also -requires implementation of the constraint functions in the model as methods called -"opening_indicator" and "sliding_indicator", see model_setup.ContactIndicators. +The solver can be further customized by specifying parameters in the model parameters. +Using the tailored line search also requires implementation of the constraint functions +in the model as methods called "opening_indicator" and "sliding_indicator", see +model_setup.ContactIndicators. """ @@ -35,17 +42,17 @@ class LineSearchNewtonSolver(pp.NewtonSolver): """Class for relaxing a nonlinear iteration update using a line search. - This class extends the iteration method to include a line search and implements - a line search based on the residual norm. + This class extends the iteration method to include a line search and implements a + line search based on the residual norm. """ def iteration(self, model) -> np.ndarray: """A single nonlinear iteration. - Add line search to the iteration method. First, call the super method to - compute a nonlinear iteration update. Then, perform a line search along - that update direction. + Add line search to the iteration method. First, call the super method to compute + a nonlinear iteration update. Then, perform a line search along that update + direction. Parameters: model: The simulation model. @@ -71,13 +78,15 @@ def nonlinear_line_search(self, model, dx: np.ndarray) -> np.ndarray: Returns: The step length vector, one for each degree of freedom. The relaxed update - is obtained by multiplying the nonlinear update step by the relaxation vector. + is obtained by multiplying the nonlinear update step by the relaxation + vector. """ return self.residual_line_search(model, dx) def residual_line_search(self, model, dx: np.ndarray) -> np.ndarray: - """Compute the relaxation factors for the current iteration based on the residual. + """Compute the relaxation factors for the current iteration based on the + residual. Parameters: model: The model. @@ -91,15 +100,14 @@ def residual_line_search(self, model, dx: np.ndarray) -> np.ndarray: return np.ones_like(dx) def objective_function(weight): - """Objective function for the trial update relaxed by - the current weight. + """Objective function for the trial update relaxed by the current weight. Parameters: weight: The relaxation factor. Returns: - The objective function value. Specified as the norm of the residual - in this implementation, may be overridden. + The objective function value. Specified as the norm of the residual in + this implementation, may be overridden. """ return self.residual_objective_function(model, dx, weight) @@ -108,9 +116,9 @@ def objective_function(weight): f_0 = objective_function(0) f_1 = objective_function(1) if f_1 < model.params["nl_convergence_tol_res"] or (f_1 < f_0 / 1e4): - # The objective function is zero at the full nonlinear step. This - # means that the nonlinear step is a minimum of the objective - # function. We can use the update without any relaxation. + # The objective function is zero at the full nonlinear step. This means that + # the nonlinear step is a minimum of the objective function. We can use the + # update without any relaxation. return np.ones_like(dx) def f_terminate(vals): @@ -143,15 +151,16 @@ def recursive_weight_from_sampling( Parameters: a, b: The interval. - condition_function: The condition function. It takes a sequence of line search - evaluations as argument and returns True if the condition is satisfied, - indicating that the recursion should be terminated. - function: The function to be tested. Returns a scalar or vector, must be compatible - with the condition function's parameter. + condition_function: The condition function. It takes a sequence of line + search evaluations as argument and returns True if the condition is + satisfied, indicating that the recursion should be terminated. + function: The function to be tested. Returns a scalar or vector, must be + compatible with the condition function's parameter. f_a: The value of the function at a. If not given, it is computed. + f_b: The value of the function at b. If not given, it is computed. num_steps: The number of sampling points in the interval [a, b]. - step_size_tolerance: The tolerance for the step size. If the step size is smaller - than this, the recursion is terminated. + step_size_tolerance: The tolerance for the step size. If the step size is + smaller than this, the recursion is terminated. Returns: The smallest tested value not satisfying the condition. @@ -174,7 +183,7 @@ def recursive_weight_from_sampling( f_a = f_c a = c else: - # There is a local minimum in the narrowed-down interval [a, c] + # There is a local minimum in the narrowed-down interval [a, c]. if step_size > step_size_tolerance: # Find it to better precision return self.recursive_weight_from_sampling( @@ -187,12 +196,12 @@ def recursive_weight_from_sampling( step_size_tolerance=step_size_tolerance, ) else: - # We're happy with the precision, return the minimum + # We're happy with the precision, return the minimum. return c - # We went through the whole interval without finding a local minimum. - # Thus, we assume that the minimum lies in [c, b]. If we have reached - # the tolerance, we return b. Otherwise, we search in [c, b]. + # We went through the whole interval without finding a local minimum. Thus, we + # assume that the minimum lies in [c, b]. If we have reached the tolerance, we + # return b. Otherwise, we search in [c, b]. if step_size < step_size_tolerance: return b else: @@ -232,9 +241,9 @@ def residual_objective_function( class SplineInterpolationLineSearch: """Class for computing relaxation factors based on spline interpolation. - This class could be seen as a tool for the technical step of performing a - line search. It also specifies that this choice should be used for the - constraint weights. + This class could be seen as a tool for the technical step of performing a line + search. It also specifies that this choice should be used for the constraint + weights. """ @@ -248,7 +257,7 @@ def compute_constraint_violation_weights( max_weight: float = 1.0, interval_target_size=1e-3, ) -> float: - """Specify that the constraint weights are computed based on spline interpolation. + """Specify that the constraint weights are computed using spline interpolation. Parameters: model: The model. @@ -258,8 +267,8 @@ def compute_constraint_violation_weights( crossing_inds: The indices where the constraint function has changed sign. f_0: The value of the constraint function at the initial point. max_weight: The maximum value for the constraint weights. - interval_target_size: The target size of the interval for the root finding. The - recursion is terminated when the interval is smaller than this. + interval_target_size: The target size of the interval for the root finding. + The recursion is terminated when the interval is smaller than this. Returns: The step length weight. @@ -267,10 +276,9 @@ def compute_constraint_violation_weights( """ if not np.any(crossing_inds): return 1.0 - # If the indicator has changed, we need to compute the relaxation - # factors. We do this by recursively narrowing down the interval until - # the interval is smaller than the tolerance using a spline - # interpolation. + # If the indicator has changed, we need to compute the relaxation factors. We do + # this by recursively narrowing down the interval until the interval is smaller + # than the tolerance using a spline interpolation. a, b = 0.0, max_weight x_0 = model.equation_system.get_variable_values(iterate_index=0) f_0 = f_0[crossing_inds] @@ -311,8 +319,8 @@ def recursive_spline_interpolation( ) -> tuple[float, float, float]: """Recursive function for finding a weight satisfying a condition. - Returns both the optimum/root (see method) and the minimal interval in - which it is assumed to lie. + Returns both the optimum/root (see method) and the minimal interval in which it + is assumed to lie. Parameters: a, b: The interval. @@ -321,16 +329,17 @@ def recursive_spline_interpolation( f_a: The value of the function at a. If not given, it is computed. f_b: The value of the function at b. If not given, it is computed. num_pts: The number of sampling points in the interval [a, b]. - step_size_tolerance: The tolerance for the step size. If the step size is smaller - than this, the recursion is terminated. - method: The method for finding the minimum of the spline. Either "minimize_scalar" - or "roots". + step_size_tolerance: The tolerance for the step size. If the step size is + smaller than this, the recursion is terminated. + method: The method for finding the minimum of the spline. Either + "minimize_scalar" or "roots". Returns: Tuple containing: The minimum of the function in the interval [a, b]. The lower bound of the interval in which the minimum is assumed to lie. The upper bound of the interval in which the minimum is assumed to lie. + """ counter = 0 while b - a > interval_target_size or counter < 1: @@ -344,7 +353,7 @@ def recursive_spline_interpolation( method=method, ) x = np.linspace(a, b, num_pts) - # Find the indices on either side of alpha + # Find the indices on either side of alpha. ind = np.searchsorted(x, alpha) if ind == 0: b = x[1] @@ -379,8 +388,8 @@ def optimum_from_spline( f_a: The value of the function at a. If not given, it is computed. f_b: The value of the function at b. If not given, it is computed. num_pts: The number of sampling points in the interval [a, b]. - method: The method for finding the minimum of the spline. Either "minimize_scalar" - or "roots". + method: The method for finding the minimum of the spline. Either + "minimize_scalar" or "roots". Returns: Tuple containing: @@ -389,7 +398,6 @@ def optimum_from_spline( The y values of the spline interpolation. """ - x = np.linspace(a, b, num_pts) y_list = [] @@ -401,10 +409,10 @@ def optimum_from_spline( else: f_pt = f(pt) if np.any(np.isnan(f_pt)): - # If we get overflow, truncate the x vector + # If we get overflow, truncate the x vector. x = x[: np.where(x == pt)[0][0]] break - # Collect function values, scalar or vector + # Collect function values, scalar or vector. y_list.append(f_pt) if isinstance(y_list[0], np.ndarray): y = np.vstack(y_list) @@ -412,6 +420,16 @@ def optimum_from_spline( y = np.array(y_list) def compute_and_postprocess_single(poly, a: float, b: float) -> float: + """Compute the minimum of one spline and postprocess the result. + + Parameters: + poly: The spline. + a, b: The interval. + + Returns: + The minimum of the spline in the interval [a, b]. + + """ if method == "minimize_scalar": minimum = scipy.optimize.minimize_scalar( lambda s: poly(s), bounds=[a, b], method="bounded" @@ -422,15 +440,14 @@ def compute_and_postprocess_single(poly, a: float, b: float) -> float: if min_x.size == 0: return b else: - # Find smallest root inside [a, b] + # Find smallest root inside [a, b]. min_x = min_x[(min_x >= a) & (min_x <= b)] if min_x.size == 0: return b else: return np.min(min_x) - # Find minima of the spline - + # Find minima of the splines. if isinstance(y_list[0], np.ndarray): all_minima = [] for i in range(y.shape[1]): @@ -446,12 +463,11 @@ def compute_and_postprocess_single(poly, a: float, b: float) -> float: class ConstraintLineSearch: - """Class for computing relaxation weights based on constraint functions - for contact mechanics. + """Class for line search based on constraint functions for contact mechanics. - The contract with the Model class is that the constraint functions are - defined in the model as Operator returning methods called "opening_indicator" - and "sliding_indicator". + The contract with the Model class is that the constraint functions are defined in + the model as Operator returning methods called "opening_indicator" and + "sliding_indicator". """ @@ -464,12 +480,8 @@ class ConstraintLineSearch: """ residual_line_search: Callable[[Any, np.ndarray], np.ndarray] - """Method for computing the relaxation factors for the current iteration - based on the residual.""" - - @property - def use_fracture_minimum(self): - return True + """Method for computing the relaxation factors for the current iteration based on + the residual.""" @property def min_line_search_weight(self): @@ -501,9 +513,9 @@ def constraint_line_search( ) -> np.ndarray: """Perform line search along the Newton step for the constraints. - This method defines the constraint weights for the contact mechanics and - how they are combined to a global weight. For more advanced combinations, - this method can be overridden or the stored weights can be accessed elsewhere. + This method defines the constraint weights for the contact mechanics and how + they are combined to a global weight. For more advanced combinations, this + method can be overridden or the stored weights can be accessed elsewhere. Parameters: model: The model. @@ -521,8 +533,8 @@ def constraint_line_search( global_weight = max_weight for sd in subdomains: sd_list = [sd] - # Compute the relaxation factors for the normal and tangential - # component of the contact mechanics. + # Compute the relaxation factors for the normal and tangential component of + # the contact mechanics. normal_weights = self.constraint_weights( model, dx, @@ -535,12 +547,12 @@ def constraint_line_search( model.sliding_indicator(sd_list), max_weight=np.minimum(max_weight, normal_weights).min(), ) - # For each cell, take minimum of tangential and normal weights + # For each cell, take minimum of tangential and normal weights. combined_weights = np.vstack((tangential_weights, normal_weights)) min_weights = np.min(combined_weights, axis=0) # Store the weights for the subdomain. This facilitates more advanced - # combinations of the constraint weights, e.g. using local weights for - # each cell. + # combinations of the constraint weights, e.g. using local weights for each + # cell. sd_weights[sd] = min_weights model.mdg.subdomain_data(sd).update({"constraint_weights": min_weights}) # Combine the weights for all subdomains to a global minimum weight. @@ -575,15 +587,15 @@ def constraint_weights( The step length vector, one for each degree of freedom of the constraint. """ - # If the sign of the function defining the regions has not changed, we - # use unitary relaxation factors + # If the sign of the function defining the regions has not changed, we use + # unitary relaxation factors. x_0 = model.equation_system.get_variable_values(iterate_index=0) violation_tol = model.params.get("constraint_violation_tolerance", 3e-1) relative_cell_tol = model.params.get( "relative_constraint_transition_tolerance", 2e-1 ) - # Compute the constraint function at the maximum weights. Specify that return value - # is an array to avoid type errors. + # Compute the constraint function at the maximum weights. Specify that return + # value is an array to avoid type errors. f_1 = cast( np.ndarray, constraint_function.value( @@ -595,7 +607,7 @@ def constraint_weights( f_0 = constraint_function.value(model.equation_system, x_0) active_inds = np.ones(f_1.shape, dtype=bool) for i in range(10): - # Only consider dofs where the constraint violation has changed sign + # Only consider dofs where the constraint violation has changed sign. violation = violation_tol * np.sign(f_1) f = constraint_function - pp.wrap_as_dense_ad_array(violation) # Absolute tolerance should be safe, as constraints are assumed to be @@ -604,8 +616,8 @@ def constraint_weights( inds = np.logical_and(np.abs(f_1) > violation_tol, f_0 * f_1 < -roundoff) if i > 0: # Ensure at least one iteration. if sum(active_inds) < max(1, relative_cell_tol * active_inds.size): - # Only a few indices are active, and the set of active indices - # does not contain any new indices. We can stop the iteration. + # Only a few indices are active, and the set of active indices does + # not contain any new indices. We can stop the iteration. break else: @@ -616,8 +628,8 @@ def constraint_weights( f_0_v = f_0 - violation # Compute the relaxation factors for the indices where the constraint - # violation has changed sign. The weight corresponds to the point at - # which the constraint function crosses zero. + # violation has changed sign. The weight corresponds to the point at which + # the constraint function crosses zero. crossing_weight = self.compute_constraint_violation_weights( model, solution_update, @@ -632,9 +644,7 @@ def constraint_weights( ) weights[inds] = weight - if not self.use_fracture_minimum: - break # Experimental. - # Check how many indices are active for current weight + # Check how many indices are active for current weight. f_1 = cast( np.ndarray, constraint_function.value( @@ -653,11 +663,3 @@ def constraint_weights( violation_tol = violation_tol / 2 return weights - - -# class NonlinearSolver( -# ConstraintLineSearch, -# SplineInterpolationLineSearch, -# LineSearchNewtonSolver, -# ): -# pass From aa38a9488cdeb61b767c5dc30758ee6cdde6aa59 Mon Sep 17 00:00:00 2001 From: Ivar Stefansson Date: Fri, 9 Aug 2024 09:50:34 +0000 Subject: [PATCH 8/8] STY: flake8 --- src/porepy/models/solution_strategy.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/porepy/models/solution_strategy.py b/src/porepy/models/solution_strategy.py index b426c05d59..7379ea39ed 100644 --- a/src/porepy/models/solution_strategy.py +++ b/src/porepy/models/solution_strategy.py @@ -853,8 +853,9 @@ def sliding_indicator( 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. + # 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)