From e1aa9f8b2ee9eba985a70f6b612d22f0d9f3e33e Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Fri, 10 Jan 2025 10:48:12 +0000 Subject: [PATCH 1/2] Transpose interpolation -> adjoint interpolation --- firedrake/__future__.py | 4 +- firedrake/assemble.py | 2 +- firedrake/interpolation.py | 108 ++++++++++-------- .../firedrake/regression/test_interp_dual.py | 15 +++ .../regression/test_interpolate_cross_mesh.py | 30 ++--- .../vertexonly/test_vertex_only_fs.py | 44 +++---- 6 files changed, 118 insertions(+), 85 deletions(-) diff --git a/firedrake/__future__.py b/firedrake/__future__.py index 713ffae96f..b0c67effdb 100644 --- a/firedrake/__future__.py +++ b/firedrake/__future__.py @@ -39,8 +39,8 @@ class CrossMeshInterpolator(Interpolator, CrossMeshInterpolatorOld): def interpolate(expr, V, *args, **kwargs): default_missing_val = kwargs.pop("default_missing_val", None) if isinstance(V, Cofunction): - transpose = bool(extract_arguments(expr)) + adjoint = bool(extract_arguments(expr)) return Interpolator( expr, V.function_space().dual(), *args, **kwargs - ).interpolate(V, transpose=transpose, default_missing_val=default_missing_val) + ).interpolate(V, adjoint=adjoint, default_missing_val=default_missing_val) return Interpolator(expr, V, *args, **kwargs).interpolate(default_missing_val=default_missing_val) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index f3049ae01c..ff9c72f4d7 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -542,7 +542,7 @@ def base_form_assembly_visitor(self, expr, tensor, *args): # Assembling the action of the Jacobian adjoint. if is_adjoint: output = tensor or firedrake.Cofunction(arg_expression[0].function_space().dual()) - return interpolator._interpolate(v, output=output, transpose=True, default_missing_val=default_missing_val) + return interpolator._interpolate(v, output=output, adjoint=True, default_missing_val=default_missing_val) # Assembling the Jacobian action. if interpolator.nargs: return interpolator._interpolate(expression, output=tensor, default_missing_val=default_missing_val) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index f27024ff28..63d45dbe6b 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -3,6 +3,7 @@ import os import tempfile import abc +import warnings import FIAT import ufl @@ -23,7 +24,7 @@ import firedrake from firedrake import tsfc_interface, utils, functionspaceimpl -from firedrake.ufl_expr import Argument, action, adjoint +from firedrake.ufl_expr import Argument, action, adjoint as expr_adjoint from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError from firedrake.petsc import PETSc from firedrake.halo import _get_mtype as get_dat_mpi_type @@ -75,7 +76,7 @@ def __init__(self, expr, v, (a) unchanged if some ``output`` is given to the :meth:`interpolate` method or (b) set to zero. Can be overwritten with the ``default_missing_val`` kwarg of :meth:`interpolate`. - This does not affect transpose interpolation. Ignored if interpolating within + This does not affect adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh` (the behaviour of a :func:`.VertexOnlyMesh` in this scenario is, at present, set when it is created). default_missing_val : float @@ -139,7 +140,7 @@ def _ufl_expr_reconstruct_(self, expr, v=None, **interp_data): # - v is a function in V (NOT an expression). # - w is a function in W. # - Maths: v = Bw -# - v_star = B.interpolate(w_star, transpose = True) +# - v_star = B.interpolate(w_star, adjoint=True) # - w_star is a cofunction in W^* (such as an assembled 1-form). # - v_star is a cofunction in V^*. # - Maths: v^* = B^* w^* @@ -175,7 +176,7 @@ def interpolate( ``True`` the corresponding values are either (a) unchanged if some ``output`` is given to the :meth:`interpolate` method or (b) set to zero. In either case, if ``default_missing_val`` is specified, that - value is used. This does not affect transpose interpolation. Ignored if + value is used. This does not affect adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh` (the behaviour of a :func:`.VertexOnlyMesh` in this scenario is, at present, set when it is created). @@ -242,7 +243,7 @@ class Interpolator(abc.ABC): ``True`` the corresponding values are either (a) unchanged if some ``output`` is given to the :meth:`interpolate` method or (b) set to zero. Can be overwritten with the ``default_missing_val`` kwarg - of :meth:`interpolate`. This does not affect transpose interpolation. + of :meth:`interpolate`. This does not affect adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh` (the behaviour of a :func:`.VertexOnlyMesh` in this scenario is, at present, set when it is created). @@ -298,7 +299,7 @@ def __init__( part=v.part())}) self.expr_renumbered = expr - def _interpolate_future(self, *function, transpose=False, default_missing_val=None): + def _interpolate_future(self, *function, transpose=None, adjoint=None, default_missing_val=None): """Define the :class:`Interpolate` object corresponding to the interpolation operation of interest. Parameters @@ -306,16 +307,18 @@ def _interpolate_future(self, *function, transpose=False, default_missing_val=No *function: firedrake.function.Function or firedrake.cofunction.Cofunction If the expression being interpolated contains an argument, then the function value to interpolate. - transpose: bool - Set to true to apply the transpose (adjoint) of the - interpolation operator. + transpose : bool + Deprecated, use adjoint instead. + adjoint: bool + Set to true to apply the adjoint of the interpolation + operator. default_missing_val: bool For interpolation across meshes: the optional value to assign to DoFs in the target mesh that are outside the source mesh. If this is not set then the values are either (a) unchanged if some ``output`` is specified to the :meth:`interpolate` method or (b) set to zero. This does not affect - transpose interpolation. Ignored if interpolating within the same + adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh`. Returns @@ -339,7 +342,10 @@ def _interpolate_future(self, *function, transpose=False, default_missing_val=No allow_missing_dofs=self._allow_missing_dofs, default_missing_val=default_missing_val) if transpose: - interp = adjoint(interp) + warnings.warn("'transpose' argument is deprecated, use 'adjoint' instead", FutureWarning) + adjoint = transpose or adjoint + if adjoint: + interp = expr_adjoint(interp) if function: f, = function @@ -349,7 +355,7 @@ def _interpolate_future(self, *function, transpose=False, default_missing_val=No return interp @PETSc.Log.EventDecorator() - def interpolate(self, *function, output=None, transpose=False, default_missing_val=None, + def interpolate(self, *function, output=None, transpose=None, adjoint=False, default_missing_val=None, ad_block_tag=None): """Compute the interpolation by assembling the appropriate :class:`Interpolate` object. @@ -360,16 +366,18 @@ def interpolate(self, *function, output=None, transpose=False, default_missing_v then the function value to interpolate. output: firedrake.function.Function or firedrake.cofunction.Cofunction A function to contain the output. - transpose: bool - Set to true to apply the transpose (adjoint) of the - interpolation operator. + transpose : bool + Deprecated, use adjoint instead. + adjoint: bool + Set to true to apply the adjoint of the interpolation + operator. default_missing_val: bool For interpolation across meshes: the optional value to assign to DoFs in the target mesh that are outside the source mesh. If this is not set then the values are either (a) unchanged if some ``output`` is specified to the :meth:`interpolate` method or (b) set to zero. This does not affect - transpose interpolation. Ignored if interpolating within the same + adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh`. ad_block_tag: str An optional string for tagging the resulting assemble block on the Pyadjoint tape. @@ -379,7 +387,6 @@ def interpolate(self, *function, output=None, transpose=False, default_missing_v firedrake.function.Function or firedrake.cofunction.Cofunction The resulting interpolated function. """ - import warnings from firedrake.assemble import assemble warnings.warn("""The use of `interpolate` to perform the numerical interpolation is deprecated. @@ -402,9 +409,12 @@ def interpolate(self, *function, output=None, transpose=False, default_missing_v Alternatively, you can also perform other symbolic operations on the interpolation operator, such as taking the derivative, and then assemble the resulting form. """, FutureWarning) + if transpose: + warnings.warn("'transpose' argument is deprecated, use 'adjoint' instead", FutureWarning) + adjoint = transpose or adjoint # Get the Interpolate object - interp = self._interpolate_future(*function, transpose=transpose, + interp = self._interpolate_future(*function, adjoint=adjoint, default_missing_val=default_missing_val) if isinstance(self.V, firedrake.Function) and not output: @@ -643,7 +653,8 @@ def _interpolate( self, *function, output=None, - transpose=False, + transpose=None, + adjoint=False, default_missing_val=None, **kwargs, ): @@ -653,9 +664,12 @@ def _interpolate( """ from firedrake.assemble import assemble - if transpose and not self.nargs: + if transpose is not None: + warnings.warn("'transpose' argument is deprecated, use 'adjoint' instead", FutureWarning) + adjoint = transpose or adjoint + if adjoint and not self.nargs: raise ValueError( - "Can currently only apply transpose interpolation with arguments." + "Can currently only apply adjoint interpolation with arguments." ) if self.nargs != len(function): raise ValueError( @@ -672,7 +686,7 @@ def _interpolate( else: f_src = self.expr - if transpose: + if adjoint: try: V_dest = self.expr.function_space().dual() except AttributeError: @@ -684,7 +698,7 @@ def _interpolate( V_dest = coeffs[0].function_space().dual() else: raise ValueError( - "Can't transpose interpolate an expression with no coefficients or arguments." + "Can't adjoint interpolate an expression with no coefficients or arguments." ) else: if isinstance(self.V, (firedrake.Function, firedrake.Cofunction)): @@ -710,14 +724,14 @@ def _interpolate( # so the sub_interpolators are already prepared to interpolate # without needing to be given a Function assert not self.nargs - interp = sub_interpolator._interpolate_future(transpose=transpose, **kwargs) + interp = sub_interpolator._interpolate_future(adjoint=adjoint, **kwargs) assemble(interp, tensor=output_sub_func) else: - interp = sub_interpolator._interpolate_future(transpose=transpose, **kwargs) + interp = sub_interpolator._interpolate_future(adjoint=adjoint, **kwargs) assemble(action(interp, f_src_sub_func), tensor=output_sub_func) return output - if not transpose: + if not adjoint: if f_src is self.expr: # f_src is already contained in self.point_eval_interpolate assert not self.nargs @@ -764,7 +778,7 @@ def _interpolate( ] = f_src_at_dest_node_coords_dest_mesh_decomp.dat.data_ro[:] else: - # adjoint/transpose interpolation + # adjoint interpolation # f_src is a cofunction on V_dest.dual as originally specified when # creating the interpolator. Our first adjoint operation is to @@ -780,25 +794,25 @@ def _interpolate( : ] = f_src.dat.data_ro[:] - # The rest of the transpose interpolation is merely the composition - # of the transpose interpolators in the reverse direction. NOTE: I + # The rest of the adjoint interpolation is merely the composition + # of the adjoint interpolators in the reverse direction. NOTE: I # don't have to worry about skipping over missing points here # because I'm going from the input ordering VOM to the original VOM # and all points from the input ordering VOM are in the original. - interp = action(adjoint(self.to_input_ordering_interpolate), f_src_at_dest_node_coords_dest_mesh_decomp) + interp = action(expr_adjoint(self.to_input_ordering_interpolate), f_src_at_dest_node_coords_dest_mesh_decomp) f_src_at_src_node_coords = assemble(interp) # NOTE: if I wanted the default missing value to be applied to - # transpose interpolation I would have to do it here. However, + # adjoint interpolation I would have to do it here. However, # this would require me to implement default missing values for - # transpose interpolation from a point evaluation interpolator + # adjoint interpolation from a point evaluation interpolator # which I haven't done. I wonder if it is necessary - perhaps the # adjoint operator always sets all the values of the resulting # cofunction? My initial attempt to insert setting the dat values - # prior to performing the multTranspose operation in + # prior to performing the multHermitian operation in # SameMeshInterpolator.interpolate did not effect the result. For # now, I say in the docstring that it only applies to forward # interpolation. - interp = action(adjoint(self.point_eval_interpolate), f_src_at_src_node_coords) + interp = action(expr_adjoint(self.point_eval_interpolate), f_src_at_src_node_coords) assemble(interp, tensor=output) return output @@ -823,13 +837,17 @@ def __init__(self, expr, V, subset=None, freeze_expr=False, access=op2.WRITE, bc self.nargs = len(arguments) @PETSc.Log.EventDecorator() - def _interpolate(self, *function, output=None, transpose=False, **kwargs): + def _interpolate(self, *function, output=None, transpose=None, adjoint=False, **kwargs): """Compute the interpolation. For arguments, see :class:`.Interpolator`. """ - if transpose and not self.nargs: - raise ValueError("Can currently only apply transpose interpolation with arguments.") + + if transpose is not None: + warnings.warn("'transpose' argument is deprecated, use 'adjoint' instead", FutureWarning) + adjoint = transpose or adjoint + if adjoint and not self.nargs: + raise ValueError("Can currently only apply adjoint interpolation with arguments.") if self.nargs != len(function): raise ValueError("Passed %d Functions to interpolate, expected %d" % (len(function), self.nargs)) @@ -851,8 +869,8 @@ def _interpolate(self, *function, output=None, transpose=False, **kwargs): function, = function if not hasattr(function, "dat"): raise ValueError("The expression had arguments: we therefore need to be given a Function (not an expression) to interpolate!") - if transpose: - mul = assembled_interpolator.handle.multTranspose + if adjoint: + mul = assembled_interpolator.handle.multHermitian V = self.arguments[0].function_space() else: mul = assembled_interpolator.handle.mult @@ -1520,7 +1538,7 @@ class VomOntoVomDummyMat(object): forward_reduce : bool If ``True``, the action of the operator (accessed via the `mult` method) is to perform a SF reduce from the source vec to the target - vec, whilst the adjoint action (accessed via the `multTranspose` + vec, whilst the adjoint action (accessed via the `multHermitian` method) is to perform a SF broadcast from the source vec to the target vec. If ``False``, the opposite is true. V : `.FunctionSpace` @@ -1630,23 +1648,23 @@ def mult(self, source_vec, target_vec): else: self.broadcast(coeff_vec, target_vec) - def multTranspose(self, source_vec, target_vec): - # can only do transpose if our expression exclusively contains a + def multHermitian(self, source_vec, target_vec): + # can only do adjoint if our expression exclusively contains a # single argument, making the application of the adjoint operator # straightforward (haven't worked out how to do this otherwise!) if not len(self.arguments) == 1: raise NotImplementedError( - "Can only apply transpose to expressions with one argument!" + "Can only apply adjoint to expressions with one argument!" ) if self.arguments[0] is not self.expr: raise NotImplementedError( - "Can only apply transpose to expressions consisting of a single argument at the moment." + "Can only apply adjoint to expressions consisting of a single argument at the moment." ) if self.forward_reduce: self.broadcast(source_vec, target_vec) else: # We need to ensure the target vec is zeroed for SF Reduce to - # represent multTranspose in case the interpolation matrix is not + # represent multHermitian in case the interpolation matrix is not # square (in which case it will have columns which are zero). This # happens when we interpolate from an input-ordering vertex-only # mesh to an immersed vertex-only mesh where the input ordering diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 5928d81e4e..404f135990 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -2,6 +2,7 @@ import numpy as np from firedrake import * from firedrake.__future__ import * +from firedrake.utils import complex_mode import ufl @@ -100,6 +101,20 @@ def test_assemble_interp_adjoint_model(V1, V2): assert np.allclose(res.dat.data, Ivfstar.dat.data) +def test_assemble_interp_adjoint_complex(mesh, V1, f1): + if complex_mode: + f1 = Constant(3 - 5.j) * f1 + + a = assemble(conj(TestFunction(V1)) * dx) + b = assemble(action(adjoint(Interpolate(f1 * TestFunction(V1), V1)), a)) + + x, y = SpatialCoordinate(mesh) + f2 = Function(V1, name="f2").interpolate( + exp(x) * y) + + assert np.allclose(assemble(b(f2)), assemble(Function(V1).interpolate(conj(f1 * f2)) * dx)) + + def test_assemble_interp_rank0(V1, V2, f1): # -- Interpolate(f1, u2) (rank 0) -- # v2 = TestFunction(V2) diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index ef9fa28769..a5f0cff6be 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -355,7 +355,7 @@ def test_interpolate_unitsquare_mixed(): assert np.allclose(got[:, 0], expected_1) assert np.allclose(got[:, 1], expected_2) cofunc_dest = assemble(inner(f_dest, TestFunction(V_dest)) * dx) - cofunc_src = assemble(interpolator.interpolate(cofunc_dest, transpose=True)) + cofunc_src = assemble(interpolator.interpolate(cofunc_dest, adjoint=True)) assert not np.allclose(f_src.dat.data_ro[0], cofunc_src.dat.data_ro[0]) assert not np.allclose(f_src.dat.data_ro[1], cofunc_src.dat.data_ro[1]) @@ -369,7 +369,7 @@ def test_interpolate_unitsquare_mixed(): def test_exact_refinement(): # With an exact mesh refinement, we can do error checks to see if our - # forward and transpose interpolations are exact where we expect them to + # forward and adjoint interpolations are exact where we expect them to # be. m_coarse = UnitSquareMesh(2, 2) m_fine = UnitSquareMesh(8, 8) @@ -397,12 +397,12 @@ def test_exact_refinement(): f_coarse_on_fine = assemble(interpolator_coarse_to_fine.interpolate(f_coarse)) assert np.allclose(f_coarse_on_fine.dat.data_ro, f_fine.dat.data_ro) - # Transpose interpolation takes us from V_fine^* to V_coarse^* so we should + # Adjoint interpolation takes us from V_fine^* to V_coarse^* so we should # also get an exact result here. cofunction_coarse = assemble(inner(f_coarse, TestFunction(V_coarse)) * dx) cofunction_fine = assemble(inner(f_fine, TestFunction(V_fine)) * dx) cofunction_fine_on_coarse = assemble(interpolator_coarse_to_fine.interpolate( - cofunction_fine, transpose=True + cofunction_fine, adjoint=True )) assert np.allclose( cofunction_fine_on_coarse.dat.data_ro, cofunction_coarse.dat.data_ro @@ -425,14 +425,14 @@ def test_exact_refinement(): f_fine_on_coarse = assemble(interpolator_fine_to_coarse.interpolate(f_fine)) assert np.allclose(f_fine_on_coarse.dat.data_ro, f_coarse.dat.data_ro) - # But transpose interpolation, which takes us from V_coarse^* to V_fine^* + # But adjoint interpolation, which takes us from V_coarse^* to V_fine^* # won't exactly reproduce the cofunction of f_coarse, since V_fine^* is a # bigger space than V_coarse^*. This only worked before because we could # exactly represent the expression in V_coarse. cofunction_fine = assemble(inner(f_fine, TestFunction(V_fine)) * dx) cofunction_coarse = assemble(inner(f_coarse, TestFunction(V_coarse)) * dx) cofunction_coarse_on_fine = assemble(interpolator_fine_to_coarse.interpolate( - cofunction_coarse, transpose=True + cofunction_coarse, adjoint=True )) assert not np.allclose( cofunction_coarse_on_fine.dat.data_ro, cofunction_fine.dat.data_ro @@ -445,10 +445,10 @@ def test_exact_refinement(): f_course_on_fine = assemble(interpolator_coarse_to_fine.interpolate(f_coarse)) assert not np.allclose(f_course_on_fine.dat.data_ro, f_fine.dat.data_ro) - # But the transpose operation, which takes us from V_fine^* to + # But the adjoint operation, which takes us from V_fine^* to # V_coarse^* correctly reproduces cofunction_coarse from cofunction_fine cofunction_fine_on_coarse = assemble(interpolator_coarse_to_fine.interpolate( - cofunction_fine, transpose=True + cofunction_fine, adjoint=True )) assert not np.allclose( cofunction_fine_on_coarse.dat.data_ro, cofunction_coarse.dat.data_ro @@ -516,7 +516,7 @@ def get_expected_values( atol, ) cofunction_dest = assemble(inner(f_dest, TestFunction(V_dest)) * dx) - interpolator_function_transpose( + interpolator_function_adjoint( interpolator, f_src, cofunction_dest, m_src, m_dest, coords, expected, atol ) interpolator_expression( @@ -597,16 +597,16 @@ def interpolator_function( # can't interpolate expressions using an interpolator assemble(interpolator.interpolate(2 * f_src)) cofunction_dest = assemble(inner(f_dest, TestFunction(V_dest)) * dx) - assemble(interpolator.interpolate(2 * cofunction_dest, transpose=True)) + assemble(interpolator.interpolate(2 * cofunction_dest, adjoint=True)) return interpolator, f_src, f_dest, m_src -def interpolator_function_transpose( +def interpolator_function_adjoint( interpolator, f_src, cofunction_dest, m_src, m_dest, coords, expected, atol ): f_dest = assemble(interpolator.interpolate(f_src)) - cofunction_dest_on_src = assemble(interpolator.interpolate(cofunction_dest, transpose=True)) + cofunction_dest_on_src = assemble(interpolator.interpolate(cofunction_dest, adjoint=True)) assert cofunction_dest_on_src.function_space().mesh() is m_src assert np.isclose( assemble(action(cofunction_dest_on_src, f_src)), @@ -634,7 +634,7 @@ def interpolator_expression( assert np.allclose(got, 2 * expected, atol=2 * atol) cofunction_dest = assemble(inner(f_dest, TestFunction(V_dest)) * dx) cofunction_dest_on_src = assemble(interpolator.interpolate( - cofunction_dest, transpose=True + cofunction_dest, adjoint=True )) assert cofunction_dest_on_src.function_space().mesh() is m_src assert np.isclose( @@ -679,7 +679,7 @@ def test_missing_dofs(): assemble(interpolator.interpolate(f_src, default_missing_val=0.0), tensor=f_dest) assert np.allclose(f_dest.at(coords), np.array([0.25, 0.0])) - # Try the other way around so we can check transpose is unaffected + # Try the other way around so we can check adjoint is unaffected m_src = UnitSquareMesh(4, 5) m_src.coordinates.dat.data[:] *= 2 m_dest = UnitSquareMesh(2, 3) @@ -692,7 +692,7 @@ def test_missing_dofs(): cofunction_src = assemble(inner(Function(V_dest), TestFunction(V_dest)) * dx) cofunction_src.dat.data_wo[:] = 1.0 cofunction_dest = assemble(interpolator.interpolate( - cofunction_src, transpose=True, default_missing_val=2.0 + cofunction_src, adjoint=True, default_missing_val=2.0 )) assert np.all(cofunction_dest.dat.data_ro != 2.0) diff --git a/tests/firedrake/vertexonly/test_vertex_only_fs.py b/tests/firedrake/vertexonly/test_vertex_only_fs.py index a900aff568..ce071267b6 100644 --- a/tests/firedrake/vertexonly/test_vertex_only_fs.py +++ b/tests/firedrake/vertexonly/test_vertex_only_fs.py @@ -139,7 +139,7 @@ def functionspace_tests(vm): g.interpolate(h) assert np.allclose(g.dat.data_ro_with_halos, np.prod(vm.coordinates.dat.data_ro_with_halos.reshape(-1, vm.geometric_dimension()), axis=1)) # Can equivalently create interpolators and use them. NOTE the - # transpose interpolator is equivilent to the inverse here because the + # adjoint interpolator is equivilent to the inverse here because the # inner product matrix in the reisz representer is the identity. TODO: when # we introduce cofunctions, this will need to be rewritten. I_io = Interpolator(TestFunction(V), W) @@ -151,27 +151,27 @@ def functionspace_tests(vm): assert np.allclose(h2.dat.data_ro_with_halos[idxs_to_include], 2*np.prod(vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include].reshape(-1, vm.input_ordering.geometric_dimension()), axis=1)) h_star = h.riesz_representation(riesz_map="l2") - g = assemble(I_io.interpolate(h_star, transpose=True)) + g = assemble(I_io.interpolate(h_star, adjoint=True)) assert np.allclose(g.dat.data_ro_with_halos, np.prod(vm.coordinates.dat.data_ro_with_halos.reshape(-1, vm.geometric_dimension()), axis=1)) with pytest.raises(NotImplementedError): - # Can't use transpose on interpolators with expressions yet - g2 = assemble(I2_io.interpolate(h_star, transpose=True)) + # Can't use adjoint on interpolators with expressions yet + g2 = assemble(I2_io.interpolate(h_star, adjoint=True)) assert np.allclose(g2.dat.data_ro_with_halos, 2*np.prod(vm.coordinates.dat.data_ro_with_halos.reshape(-1, vm.geometric_dimension()), axis=1)) - I_io_transpose = Interpolator(TestFunction(W), V) - I2_io_transpose = Interpolator(2*TestFunction(W), V) - h_star = assemble(I_io_transpose.interpolate(g, transpose=True)) + I_io_adjoint = Interpolator(TestFunction(W), V) + I2_io_adjoint = Interpolator(2*TestFunction(W), V) + h_star = assemble(I_io_adjoint.interpolate(g, adjoint=True)) h = h_star.riesz_representation(riesz_map="l2") assert np.allclose(h.dat.data_ro_with_halos[idxs_to_include], np.prod(vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include].reshape(-1, vm.input_ordering.geometric_dimension()), axis=1)) assert np.all(h.dat.data_ro_with_halos[~idxs_to_include] == 0) with pytest.raises(NotImplementedError): - # Can't use transpose on interpolators with expressions yet - h2 = assemble(I2_io_transpose.interpolate(g, transpose=True)) + # Can't use adjoint on interpolators with expressions yet + h2 = assemble(I2_io_adjoint.interpolate(g, adjoint=True)) assert np.allclose(h2.dat.data_ro_with_halos[idxs_to_include], 2*np.prod(vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include].reshape(-1, vm.input_ordering.geometric_dimension()), axis=1)) - g = assemble(I_io_transpose.interpolate(h)) + g = assemble(I_io_adjoint.interpolate(h)) assert np.allclose(g.dat.data_ro_with_halos, np.prod(vm.coordinates.dat.data_ro_with_halos.reshape(-1, vm.geometric_dimension()), axis=1)) - g2 = assemble(I2_io_transpose.interpolate(h)) + g2 = assemble(I2_io_adjoint.interpolate(h)) assert np.allclose(g2.dat.data_ro_with_halos, 2*np.prod(vm.coordinates.dat.data_ro_with_halos.reshape(-1, vm.geometric_dimension()), axis=1)) @@ -244,7 +244,7 @@ def vectorfunctionspace_tests(vm): g.interpolate(h) assert np.allclose(g.dat.data_ro_with_halos, 2*vm.coordinates.dat.data_ro_with_halos) # Can equivalently create interpolators and use them. NOTE the - # transpose interpolator is equivilent to the inverse here because the + # adjoint interpolator is equivilent to the inverse here because the # inner product matrix in the reisz representer is the identity. TODO: when # we introduce cofunctions, this will need to be rewritten. I_io = Interpolator(TestFunction(V), W) @@ -256,27 +256,27 @@ def vectorfunctionspace_tests(vm): assert np.allclose(h2.dat.data_ro[idxs_to_include], 4*vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include]) h_star = h.riesz_representation(riesz_map="l2") - g = assemble(I_io.interpolate(h_star, transpose=True)) + g = assemble(I_io.interpolate(h_star, adjoint=True)) assert np.allclose(g.dat.data_ro_with_halos, 2*vm.coordinates.dat.data_ro_with_halos) with pytest.raises(NotImplementedError): - # Can't use transpose on interpolators with expressions yet - g2 = assemble(I2_io.interpolate(h_star, transpose=True)) + # Can't use adjoint on interpolators with expressions yet + g2 = assemble(I2_io.interpolate(h_star, adjoint=True)) assert np.allclose(g2.dat.data_ro_with_halos, 4*vm.coordinates.dat.data_ro_with_halos) - I_io_transpose = Interpolator(TestFunction(W), V) - I2_io_transpose = Interpolator(2*TestFunction(W), V) - h_star = assemble(I_io_transpose.interpolate(g, transpose=True)) + I_io_adjoint = Interpolator(TestFunction(W), V) + I2_io_adjoint = Interpolator(2*TestFunction(W), V) + h_star = assemble(I_io_adjoint.interpolate(g, adjoint=True)) assert np.allclose(h_star.dat.data_ro[idxs_to_include], 2*vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include]) assert np.all(h_star.dat.data_ro_with_halos[~idxs_to_include] == 0) with pytest.raises(NotImplementedError): - # Can't use transpose on interpolators with expressions yet - h2 = assemble(I2_io_transpose.interpolate(g, transpose=True)) + # Can't use adjoint on interpolators with expressions yet + h2 = assemble(I2_io_adjoint.interpolate(g, adjoint=True)) assert np.allclose(h2.dat.data_ro[idxs_to_include], 4*vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include]) h = h_star.riesz_representation(riesz_map="l2") - g = assemble(I_io_transpose.interpolate(h)) + g = assemble(I_io_adjoint.interpolate(h)) assert np.allclose(g.dat.data_ro_with_halos, 2*vm.coordinates.dat.data_ro_with_halos) - g2 = assemble(I2_io_transpose.interpolate(h)) + g2 = assemble(I2_io_adjoint.interpolate(h)) assert np.allclose(g2.dat.data_ro_with_halos, 4*vm.coordinates.dat.data_ro_with_halos) From c4d0a65b6bac8fd91016e8ea9b243e4e8df1dd48 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Fri, 10 Jan 2025 11:08:51 +0000 Subject: [PATCH 2/2] Fix --- firedrake/interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 63d45dbe6b..c4797011dd 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -299,7 +299,7 @@ def __init__( part=v.part())}) self.expr_renumbered = expr - def _interpolate_future(self, *function, transpose=None, adjoint=None, default_missing_val=None): + def _interpolate_future(self, *function, transpose=None, adjoint=False, default_missing_val=None): """Define the :class:`Interpolate` object corresponding to the interpolation operation of interest. Parameters