Remove unused interpolation code

Author: jrmaddison

fenics_ice/interpolation.py

@@ -2,175 +2,22 @@
 """
 
 from .backend import (
-    Cell, Mesh, MeshEditor, Point, UserExpression, backend_Constant,
-    backend_Function, backend_ScalarType, parameters)
+    Cell, Mesh, MeshEditor, Point, backend_ScalarType, parameters)
 from ..interface import (
-    check_space_type, comm_dup_cached, packed, space_comm, space_eq,
-    var_assign, var_comm, var_get_values, var_id, var_inner, var_is_scalar,
-    var_local_size, var_new, var_new_conjugate_dual, var_replacement,
-    var_scalar_value, var_set_values)
+    comm_dup_cached, space_comm, var_comm, var_get_values, var_local_size,
+    var_new, var_set_values)
 
-from ..equation import Equation, ZeroAssignment
-from ..equations import MatrixActionRHS
-from ..linear_equation import LinearEquation, Matrix
-from ..manager import manager_disabled
-
-from .expr import (
-    ExprEquation, derivative, eliminate_zeros, expr_zero, extract_dependencies,
-    extract_variables)
-from .variables import ReplacementConstant
+from ..linear_equation import Matrix
 
 import functools
 import mpi4py.MPI as MPI
 import numpy as np
-try:
-    import ufl_legacy as ufl
-except ModuleNotFoundError:
-    import ufl
 
 __all__ = \
     [
-        "ExprInterpolation",
-        "Interpolation",
-        "PointInterpolation"
     ]
 
 
-@manager_disabled()
-def interpolate_expression(x, expr, *, adj_x=None):
-    if adj_x is None:
-        check_space_type(x, "primal")
-    else:
-        check_space_type(x, "conjugate_dual")
-        check_space_type(adj_x, "conjugate_dual")
-    for dep in extract_variables(expr):
-        check_space_type(dep, "primal")
-
-    expr = eliminate_zeros(expr)
-
-    class Expr(UserExpression):
-        def eval(self, value, x):
-            value[:] = expr(tuple(x))
-
-        def value_shape(self):
-            return x.ufl_shape
-
-    if adj_x is None:
-        if isinstance(x, backend_Constant):
-            if isinstance(expr, backend_Constant):
-                value = expr
-            else:
-                if len(x.ufl_shape) > 0:
-                    raise ValueError("Scalar Constant required")
-                value = x.values()
-                Expr().eval(value, ())
-                value, = value
-            var_assign(x, value)
-        elif isinstance(x, backend_Function):
-            try:
-                x.assign(expr)
-            except RuntimeError:
-                x.interpolate(Expr())
-        else:
-            raise TypeError(f"Unexpected type: {type(x)}")
-    else:
-        expr_val = var_new_conjugate_dual(adj_x)
-        expr_arguments = ufl.algorithms.extract_arguments(expr)
-        if len(expr_arguments) > 0:
-            test, = expr_arguments
-            if len(test.ufl_shape) > 0:
-                raise NotImplementedError("Case not implemented")
-            expr = ufl.replace(expr, {test: ufl.classes.IntValue(1)})
-        interpolate_expression(expr_val, expr)
-
-        if isinstance(x, backend_Constant):
-            if len(x.ufl_shape) > 0:
-                raise ValueError("Scalar Constant required")
-            var_assign(x, var_inner(adj_x, expr_val))
-        elif isinstance(x, backend_Function):
-            if not space_eq(x.function_space(), adj_x.function_space()):
-                raise ValueError("Unable to perform transpose interpolation")
-            var_set_values(
-                x, var_get_values(expr_val).conjugate() * var_get_values(adj_x))  # noqa: E501
-        else:
-            raise TypeError(f"Unexpected type: {type(x)}")
-
-
-class ExprInterpolation(ExprEquation):
-    r"""Represents interpolation of `rhs` onto the space for `x`.
-
-    The forward residual :math:`\mathcal{F}` is defined so that :math:`\partial
-    \mathcal{F} / \partial x` is the identity.
-
-    :arg x: The forward solution.
-    :arg rhs: A :class:`ufl.core.expr.Expr` defining the expression to
-        interpolate onto the space for `x`. Should not depend on `x`.
-    """
-
-    def __init__(self, x, rhs):
-        deps, nl_deps = extract_dependencies(rhs, space_type="primal")
-        if var_id(x) in deps:
-            raise ValueError("Invalid dependency")
-        deps, nl_deps = list(deps.values()), tuple(nl_deps.values())
-        deps.insert(0, x)
-
-        super().__init__(x, deps, nl_deps=nl_deps, ic=False, adj_ic=False)
-        self._rhs = rhs
-
-    def drop_references(self):
-        replace_map = {dep: var_replacement(dep)
-                       for dep in self.dependencies()}
-
-        super().drop_references()
-        self._rhs = ufl.replace(self._rhs, replace_map)
-
-    def forward_solve(self, x, deps=None):
-        interpolate_expression(x, self._replace(self._rhs, deps))
-
-    def adjoint_derivative_action(self, nl_deps, dep_index, adj_x):
-        eq_deps = self.dependencies()
-        if dep_index <= 0 or dep_index >= len(eq_deps):
-            raise ValueError("Unexpected dep_index")
-
-        dep = eq_deps[dep_index]
-
-        if isinstance(dep, (backend_Constant, ReplacementConstant)):
-            if len(dep.ufl_shape) > 0:
-                raise NotImplementedError("Case not implemented")
-            dF = derivative(self._rhs, dep, argument=ufl.classes.IntValue(1))
-        else:
-            dF = derivative(self._rhs, dep)
-        dF = eliminate_zeros(dF)
-        dF = self._nonlinear_replace(dF, nl_deps)
-
-        F = var_new_conjugate_dual(dep)
-        interpolate_expression(F, dF, adj_x=adj_x)
-        return (-1.0, F)
-
-    def adjoint_jacobian_solve(self, adj_x, nl_deps, b):
-        return b
-
-    def tangent_linear(self, tlm_map):
-        x = self.x()
-
-        tlm_rhs = expr_zero(x)
-        for dep in self.dependencies():
-            if dep != x:
-                tau_dep = tlm_map[dep]
-                if tau_dep is not None:
-                    tlm_rhs = (tlm_rhs
-                               + derivative(self._rhs, dep, argument=tau_dep))
-
-        if isinstance(tlm_rhs, ufl.classes.Zero):
-            return ZeroAssignment(tlm_map[x])
-        else:
-            return ExprInterpolation(tlm_map[x], tlm_rhs)
-
-
-def function_coords(x):
-    return x.function_space().tabulate_dof_coordinates()
-
-
 def has_ghost_cells(mesh):
     for cell in range(mesh.num_cells()):
         if Cell(mesh, cell).is_ghost():
@@ -436,165 +283,3 @@ def adjoint_action(self, nl_deps, adj_x, b, b_index=0, *, method="assign"):
             b.vector()[:] -= self._P_T.dot(var_get_values(adj_x))
         else:
             raise ValueError(f"Invalid method: '{method:s}'")
-
-
-class Interpolation(LinearEquation):
-    r"""Represents interpolation of the scalar-valued function `y` onto the
-    space for `x`.
-
-    The forward residual :math:`\mathcal{F}` is defined so that :math:`\partial
-    \mathcal{F} / \partial x` is the identity.
-
-    Internally this builds (or uses a supplied) interpolation matrix for the
-    local process *only*. This behaves correctly if the there are no edges
-    between owned and non-owned nodes in the degree of freedom graph associated
-    with the discrete function space for `y`.
-
-    :arg x: A scalar-valued DOLFIN `Function` defining the forward solution.
-    :arg y: A scalar-valued DOLFIN `Function` to interpolate onto the space for
-        `x`.
-    :arg X_coords: A :class:`numpy.ndarray` defining the coordinates at which
-        to interpolate `y`. Shape is `(n, d)` where `n` is the number of
-        process local degrees of freedom for `x` and `d` is the geometric
-        dimension. Defaults to the process local degree of freedom locations
-        for `x`. Ignored if `P` is supplied.
-    :arg P: The interpolation matrix. A :class:`scipy.sparse.spmatrix`.
-    :arg tolerance: Maximum permitted distance (as returned by the DOLFIN
-        `BoundingBoxTree.compute_closest_entity` method) of an interpolation
-        point from a cell in the mesh for `y`. Ignored if `P` is supplied.
-    """
-
-    def __init__(self, x, y, *, x_coords=None, P=None,
-                 tolerance=0.0):
-        check_space_type(x, "primal")
-        check_space_type(y, "primal")
-
-        if not isinstance(x, backend_Function):
-            raise TypeError("Solution must be a Function")
-        if len(x.ufl_shape) > 0:
-            raise ValueError("Solution must be a scalar-valued Function")
-        if not isinstance(y, backend_Function):
-            raise TypeError("y must be a Function")
-        if len(y.ufl_shape) > 0:
-            raise ValueError("y must be a scalar-valued Function")
-        if (x_coords is not None) and (var_comm(x).size > 1):
-            raise TypeError("Cannot prescribe x_coords in parallel")
-
-        if P is None:
-            y_space = y.function_space()
-
-            if x_coords is None:
-                x_coords = function_coords(x)
-
-            y_cells, y_distances = point_cells(x_coords, y_space.mesh())
-            if (y_distances > tolerance).any():
-                raise RuntimeError("Unable to locate one or more cells")
-
-            y_colors = greedy_coloring(y_space)
-            P = interpolation_matrix(x_coords, y, y_cells, y_colors)
-        else:
-            P = P.copy()
-
-        super().__init__(
-            x, MatrixActionRHS(LocalMatrix(P), y))
-
-
-class PointInterpolation(Equation):
-    r"""Represents interpolation of a scalar-valued function at given points.
-
-    The forward residual :math:`\mathcal{F}` is defined so that :math:`\partial
-    \mathcal{F} / \partial x` is the identity.
-
-    Internally this builds (or uses a supplied) interpolation matrix for the
-    local process *only*. This behaves correctly if the there are no edges
-    between owned and non-owned nodes in the degree of freedom graph associated
-    with the discrete function space for `y`.
-
-    :arg X: A scalar variable, or a :class:`Sequence` of scalar variables,
-        defining the forward solution.
-    :arg y: A scalar-valued DOLFIN `Function` to interpolate.
-    :arg X_coords: A :class:`numpy.ndarray` defining the coordinates at which
-        to interpolate `y`. Shape is `(n, d)` where `n` is the number of
-        interpolation points and `d` is the geometric dimension. Ignored if `P`
-        is supplied.
-    :arg P: The interpolation matrix. A :class:`scipy.sparse.spmatrix`.
-    :arg tolerance: Maximum permitted distance (as returned by the DOLFIN
-        `BoundingBoxTree.compute_closest_entity` method) of an interpolation
-        point from a cell in the mesh for `y`. Ignored if `P` is supplied.
-    """
-
-    def __init__(self, X, y, X_coords=None, *,
-                 P=None, tolerance=0.0):
-        X = packed(X)
-        for x in X:
-            check_space_type(x, "primal")
-            if not var_is_scalar(x):
-                raise ValueError("Solution must be a scalar variable, or a "
-                                 "Sequence of scalar variables")
-        check_space_type(y, "primal")
-
-        if X_coords is None:
-            if P is None:
-                raise TypeError("X_coords required when P is not supplied")
-        else:
-            if len(X) != X_coords.shape[0]:
-                raise ValueError("Invalid number of variables")
-        if not isinstance(y, backend_Function):
-            raise TypeError("y must be a Function")
-        if len(y.ufl_shape) > 0:
-            raise ValueError("y must be a scalar-valued Function")
-
-        if P is None:
-            y_space = y.function_space()
-
-            y_cells = point_owners(X_coords, y_space, tolerance=tolerance)
-            y_colors = greedy_coloring(y_space)
-            P = interpolation_matrix(X_coords, y, y_cells, y_colors)
-        else:
-            P = P.copy()
-
-        super().__init__(X, list(X) + [y], nl_deps=[], ic=False, adj_ic=False)
-        self._P = P
-        self._P_T = P.T
-
-    def forward_solve(self, X, deps=None):
-        y = (self.dependencies() if deps is None else deps)[-1]
-
-        check_space_type(y, "primal")
-        y_v = var_get_values(y)
-        x_v_local = np.full(len(X), np.nan, dtype=backend_ScalarType)
-        for i in range(len(X)):
-            x_v_local[i] = self._P.getrow(i).dot(y_v)
-
-        comm = var_comm(y)
-        x_v = np.full(len(X), np.nan, dtype=x_v_local.dtype)
-        comm.Allreduce(x_v_local, x_v, op=MPI.SUM)
-
-        for i, x in enumerate(X):
-            var_assign(x, x_v[i])
-
-    def adjoint_derivative_action(self, nl_deps, dep_index, adj_X):
-        if dep_index != len(self.X()):
-            raise ValueError("Unexpected dep_index")
-
-        adj_x_v = np.full(len(adj_X), np.nan, dtype=backend_ScalarType)
-        for i, adj_x in enumerate(adj_X):
-            adj_x_v[i] = var_scalar_value(adj_x)
-
-        F = var_new_conjugate_dual(self.dependencies()[-1])
-        var_set_values(F, self._P_T.dot(adj_x_v))
-        return (-1.0, F)
-
-    def adjoint_jacobian_solve(self, adj_X, nl_deps, B):
-        return B
-
-    def tangent_linear(self, tlm_map):
-        X = self.X()
-        y = self.dependencies()[-1]
-
-        tlm_y = tlm_map[y]
-        if tlm_y is None:
-            return ZeroAssignment([tlm_map[x] for x in X])
-        else:
-            return PointInterpolation([tlm_map[x] for x in X], tlm_y,
-                                      P=self._P)

fenics_ice/solver.py

@@ -55,7 +55,7 @@ def interior(x_coords, y_space):
 
 
 def interpolation_matrix(x_coords, y_space):
-    from tlm_adjoint.fenics.interpolation import (
+    from .interpolation import (
         greedy_coloring, interpolation_matrix, point_owners)
 
     y_cells = point_owners(x_coords, y_space, tolerance=np.inf)
@@ -1292,7 +1292,7 @@ def comp_J_inv(self, verbose=False):
             cache_jacobian=False, cache_rhs_assembly=False).solve()
 
         if not hasattr(self, "_cached_J_mismatch_data"):
-            from tlm_adjoint.fenics.interpolation import LocalMatrix
+            from .interpolation import LocalMatrix
             from scipy.sparse import spdiags
 
             obs_local, P = interpolation_matrix(uv_obs_pts, interp_space)