diff --git a/animate/interpolation.py b/animate/interpolation.py index 8f9fc2f3..5f77102f 100644 --- a/animate/interpolation.py +++ b/animate/interpolation.py @@ -24,7 +24,7 @@ @PETSc.Log.EventDecorator() def transfer(source, target_space, transfer_method="project", **kwargs): r""" - Overload functions :func:`firedrake.__future__.interpolate` and + Overload functions :func:`firedrake.interpolation.interpolate` and :func:`firedrake.projection.project` to account for the case of two mixed function spaces defined on different meshes and for the adjoint interpolation operator when applied to :class:`firedrake.cofunction.Cofunction`\s. @@ -46,7 +46,7 @@ def transfer(source, target_space, transfer_method="project", **kwargs): :rtype: :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction` - Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` or + Extra keyword arguments are passed to :func:`firedrake.interpolation.interpolate` or :func:`firedrake.projection.project`. """ if transfer_method not in ("interpolate", "project"): @@ -66,63 +66,53 @@ def transfer(source, target_space, transfer_method="project", **kwargs): raise TypeError( "Second argument must be a FunctionSpace, Function, or Cofunction." ) - if isinstance(source, firedrake.Cofunction): - return _transfer_adjoint(source, target, transfer_method, **kwargs) - elif source.function_space() == target.function_space(): - return target.assign(source) + if transfer_method == "interpolate": + return interpolate(source, target, **kwargs) else: - return _transfer_forward(source, target, transfer_method, **kwargs) + return project(source, target, **kwargs) -@PETSc.Log.EventDecorator() -def interpolate(source, target_space, **kwargs): - r""" - Overload function :func:`firedrake.__future__.interpolate` to account for the case - of two mixed function spaces defined on different meshes and for the adjoint - interpolation operator when applied to :class:`firedrake.cofunction.Cofunction`\s. - - :arg source: the function to be transferred - :type source: :class:`firedrake.function.Function` or - :class:`firedrake.cofunction.Cofunction` - :arg target_space: the function space which we seek to transfer onto, or the - function or cofunction to use as the target - :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`, - :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction` - :returns: the transferred function - :rtype: :class:`firedrake.function.Function` or - :class:`firedrake.cofunction.Cofunction` - - Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` - """ - return transfer(source, target_space, transfer_method="interpolate", **kwargs) +def _validate_consistent_spaces(Vs, Vt): + if Vs._dual != Vt._dual: + raise ValueError("Spaces must be both primal or both dual.") + if hasattr(Vs, "num_sub_spaces"): + if not hasattr(Vt, "num_sub_spaces"): + raise ValueError( + "Source space has multiple components but target space does not." + ) + if Vs.num_sub_spaces() != Vt.num_sub_spaces(): + raise ValueError( + "Inconsistent numbers of components in source and target spaces:" + f" {Vs.num_sub_spaces()} vs. {Vt.num_sub_spaces()}." + ) + elif hasattr(Vt, "num_sub_spaces"): + raise ValueError( + "Target space has multiple components but source space does not." + ) @PETSc.Log.EventDecorator() -def project(source, target_space, **kwargs): +def interpolate(source, target, **kwargs): r""" - Overload function :func:`firedrake.projection.project` to account for the case of - two mixed function spaces defined on different meshes and for the adjoint - projection operator when applied to :class:`firedrake.cofunction.Cofunction`\s. - - For details on the approach for achieving boundedness through mass lumping and - post-processing, see :cite:`Farrell:2009`. + Overload :func:`firedrake.interpolation.interpolate` to account for the case of + mixed function spaces. - :arg source: the function to be transferred + :arg source: the function or cofunction to be transferred :type source: :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction` - :arg target_space: the function space which we seek to transfer onto, or the - function or cofunction to use as the target - :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`, - :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction` - :returns: the transferred function - :rtype: :class:`firedrake.function.Function` or + :arg target: the function or cofunction to use as the target, which is modified in + place + :type target: :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction` - :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness - :type bounded: :class:`bool` - Extra keyword arguments are passed to :func:`firedrake.projection.project`. + Extra keyword arguments are passed to :func:`firedrake.interpolation.interpolate` """ - return transfer(source, target_space, transfer_method="project", **kwargs) + _validate_consistent_spaces(source.function_space(), target.function_space()) + if hasattr(target.function_space(), "num_sub_spaces"): + for s, t in zip(source.subfunctions, target.subfunctions, strict=True): + t.interpolate(s, **kwargs) + else: + target.interpolate(source, **kwargs) # TODO: Reimplement by introducing a LumpedSupermeshProjector subclass of @@ -148,148 +138,46 @@ def _supermesh_project(source, target, bounded=False): @PETSc.Log.EventDecorator() -def _transfer_forward(source, target, transfer_method, **kwargs): - """ - Apply mesh-to-mesh transfer operator to a Function. +def project(source, target, bounded=False, **kwargs): + r""" + Overload :func:`firedrake.projection.project` to account for the case of mixed + function spaces and handle dual Functions, with a bounded option. - This function extends the functionality of :func:`firedrake.__future__.interpolate` - and :func:`firedrake.projection.project` to account for mixed spaces. + For details on the approach for achieving boundedness through mass lumping and + post-processing, see :cite:`Farrell:2009`. - :arg source: the Function to be transferred - :type source: :class:`firedrake.function.Function` - :arg target: the Function which we seek to transfer onto - :type target: :class:`firedrake.function.Function` - :kwarg transfer_method: the method to use for the transfer. Options are - "interpolate" (default) and "project" - :type transfer_method: str + :arg source: the function or cofunction to be transferred + :type source: :class:`firedrake.function.Function` or + :class:`firedrake.cofunction.Cofunction` + :arg target: the function or cofunction to transfer onto, which is modified in + place + :type target: :class:`firedrake.function.Function` or + :class:`firedrake.cofunction.Cofunction` :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness - (project only) :type bounded: :class:`bool` - :returns: the transferred Function - :rtype: :class:`firedrake.function.Function` - Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` or - :func:`firedrake.projection.project`. + Extra keyword arguments are passed to :func:`firedrake.projection.project`. """ - is_project = transfer_method == "project" - bounded = is_project and kwargs.pop("bounded", False) - Vs = source.function_space() - Vt = target.function_space() - _validate_matching_spaces(Vs, Vt) - assert isinstance(target, firedrake.Function) - if hasattr(Vt, "num_sub_spaces"): - for s, t in zip(source.subfunctions, target.subfunctions, strict=False): - if transfer_method == "interpolate": - t.interpolate(s, **kwargs) - elif transfer_method == "project": - if bounded: - _supermesh_project(s, t, bounded=True) - else: - t.project(s, **kwargs) - else: - raise ValueError( - f"Invalid transfer method: {transfer_method}." - " Options are 'interpolate' and 'project'." - ) - else: - if transfer_method == "interpolate": - target.interpolate(source, **kwargs) - elif transfer_method == "project": + _validate_consistent_spaces(source.function_space(), target.function_space()) + if not source.function_space()._dual: + for s, t in zip(source.subfunctions, target.subfunctions, strict=True): if bounded: - _supermesh_project(source, target, bounded=True) + _supermesh_project(s, t, bounded=True) else: - target.project(source, **kwargs) - else: - raise ValueError( - f"Invalid transfer method: {transfer_method}." - " Options are 'interpolate' and 'project'." - ) - return target - - -@PETSc.Log.EventDecorator() -def _transfer_adjoint(target_b, source_b, transfer_method, **kwargs): - """ - Apply an adjoint mesh-to-mesh transfer operator to a Cofunction. - - :arg target_b: seed Cofunction from the target space of the forward projection - :type target_b: :class:`firedrake.cofunction.Cofunction` - :arg source_b: output Cofunction from the source space of the forward projection - :type source_b: :class:`firedrake.cofunction.Cofunction` - :kwarg transfer_method: the method to use for the transfer. Options are - "interpolate" (default) and "project" - :type transfer_method: str - :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness - (project only) - :type bounded: :class:`bool` - :returns: the transferred Cofunction - :rtype: :class:`firedrake.cofunction.Cofunction` - - Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` or - :func:`firedrake.projection.project`. - """ - is_project = transfer_method == "project" - bounded = is_project and kwargs.pop("bounded", False) - - # Map to Functions to apply the adjoint transfer - if not isinstance(target_b, firedrake.Function): - target_b = cofunction2function(target_b) - if not isinstance(source_b, firedrake.Function): - source_b = cofunction2function(source_b) - - Vt = target_b.function_space() - Vs = source_b.function_space() - if Vs == Vt: - source_b.assign(target_b) - return function2cofunction(source_b) - - _validate_matching_spaces(Vs, Vt) - if hasattr(Vs, "num_sub_spaces"): - target_b_split = target_b.subfunctions - source_b_split = source_b.subfunctions + t.project(s, **kwargs) else: - target_b_split = [target_b] - source_b_split = [source_b] - - # Apply adjoint transfer operator to each component - for i, (t_b, s_b) in enumerate(zip(target_b_split, source_b_split, strict=False)): - if transfer_method == "interpolate": - raise NotImplementedError( - "Adjoint of interpolation operator not implemented." - ) # TODO (#113) - elif transfer_method == "project": + for s, t in zip(source.subfunctions, target.subfunctions, strict=True): + sf = cofunction2function(s) + tf = cofunction2function(t) + Vs = sf.function_space() ksp = petsc4py.KSP().create() - ksp.setOperators(assemble_mass_matrix(t_b.function_space(), lumped=bounded)) - mixed_mass = assemble_mixed_mass_matrix(Vt[i], Vs[i]) - with t_b.dat.vec_ro as tb, s_b.dat.vec_wo as sb: - residual = tb.copy() - ksp.solveTranspose(tb, residual) - mixed_mass.mult(residual, sb) # NOTE: already transposed above - else: - raise ValueError( - f"Invalid transfer method: {transfer_method}." - " Options are 'interpolate' and 'project'." - ) - - # Map back to a Cofunction - return function2cofunction(source_b) - - -def _validate_matching_spaces(Vs, Vt): - if hasattr(Vs, "num_sub_spaces"): - if not hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Source space has multiple components but target space does not." - ) - if Vs.num_sub_spaces() != Vt.num_sub_spaces(): - raise ValueError( - "Inconsistent numbers of components in source and target spaces:" - f" {Vs.num_sub_spaces()} vs. {Vt.num_sub_spaces()}." - ) - elif hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Target space has multiple components but source space does not." - ) + ksp.setOperators(assemble_mass_matrix(Vs, lumped=bounded)) + mixed_mass = assemble_mixed_mass_matrix(Vs, tf.function_space()) + with sf.dat.vec_ro as vs, tf.dat.vec_wo as vt: + residual = vs.copy() + ksp.solveTranspose(vs, residual) + mixed_mass.mult(residual, vt) # NOTE: already transposed above + function2cofunction(tf, cofunc=t) @PETSc.Log.EventDecorator() diff --git a/animate/metric.py b/animate/metric.py index b6702d32..2c2be202 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -7,7 +7,7 @@ import numpy as np import sympy import ufl -from firedrake.__future__ import interpolate +from firedrake.interpolation import interpolate from firedrake.petsc import PETSc from pyop2 import op2 diff --git a/animate/quality.py b/animate/quality.py index d3e060f4..fc47226c 100644 --- a/animate/quality.py +++ b/animate/quality.py @@ -6,7 +6,7 @@ import firedrake import ufl -from firedrake.__future__ import interpolate +from firedrake.interpolation import interpolate from firedrake.petsc import PETSc from pyop2 import op2 from pyop2.utils import get_petsc_dir diff --git a/animate/recovery.py b/animate/recovery.py index eb833589..051c9507 100644 --- a/animate/recovery.py +++ b/animate/recovery.py @@ -6,7 +6,7 @@ import firedrake import ufl -from firedrake.__future__ import interpolate +from firedrake.interpolation import interpolate from firedrake.petsc import PETSc from pyop2 import op2 diff --git a/animate/utility.py b/animate/utility.py index caafb7a7..0dacf50a 100644 --- a/animate/utility.py +++ b/animate/utility.py @@ -8,7 +8,7 @@ import firedrake.function as ffunc import firedrake.mesh as fmesh import ufl -from firedrake.__future__ import interpolate +from firedrake.interpolation import interpolate from firedrake.petsc import PETSc from mpi4py import MPI @@ -291,14 +291,17 @@ def assemble_mass_matrix(space, norm_type="L2", lumped=False): return mass_matrix.createDiagonal(x) -def cofunction2function(cofunc): +def cofunction2function(cofunc, func=None): """ :arg cofunc: a cofunction :type cofunc: :class:`firedrake.cofunction.Cofunction` + :kwarg func: a function for the return value + :type func: :class:`firedrake.function.Function` :returns: a function with the same underyling data :rtype: :class:`firedrake.function.Function` """ - func = ffunc.Function(cofunc.function_space().dual()) + if func is None: + func = ffunc.Function(cofunc.function_space().dual()) if isinstance(func.dat.data_with_halos, tuple): for i, arr in enumerate(func.dat.data_with_halos): arr[:] = cofunc.dat.data_with_halos[i] @@ -307,14 +310,17 @@ def cofunction2function(cofunc): return func -def function2cofunction(func): +def function2cofunction(func, cofunc=None): """ :arg func: a function :type func: :class:`firedrake.function.Function` + :kwarg cofunc: a cofunction for the return value + :type cofunc: :class:`firedrake.cofunction.Cofunction` :returns: a cofunction with the same underlying data :rtype: :class:`firedrake.cofunction.Cofunction` """ - cofunc = firedrake.Cofunction(func.function_space().dual()) + if cofunc is None: + cofunc = firedrake.Cofunction(func.function_space().dual()) if isinstance(cofunc.dat.data_with_halos, tuple): for i, arr in enumerate(cofunc.dat.data_with_halos): arr[:] = func.dat.data_with_halos[i] diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 5f22aa07..75118b4a 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -22,13 +22,11 @@ from animate.interpolation import ( _supermesh_project, - _transfer_adjoint, - _transfer_forward, clement_interpolant, project, transfer, ) -from animate.utility import function2cofunction +from animate.utility import cofunction2function, function2cofunction class TestClement(unittest.TestCase): @@ -168,15 +166,6 @@ def test_method_typo_error(self): msg = "Invalid transfer method: proj. Options are 'interpolate' and 'project'." self.assertEqual(str(cm.exception), msg) - @parameterized.expand([_transfer_forward, _transfer_adjoint]) - def test_method_typo_error_private(self, private_transfer_func): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - with self.assertRaises(ValueError) as cm: - private_transfer_func(Function(Vs), Function(Vt), transfer_method="proj") - msg = "Invalid transfer method: proj. Options are 'interpolate' and 'project'." - self.assertEqual(str(cm.exception), msg) - @parameterized.expand(["interpolate", "project"]) def test_notimplemented_error(self, transfer_method): Vs = FunctionSpace(self.source_mesh, "CG", 1) @@ -186,6 +175,13 @@ def test_notimplemented_error(self, transfer_method): msg = f"Can only currently {transfer_method} Functions and Cofunctions." self.assertEqual(str(cm.exception), msg) + def test_primal_dual_inconsistency_error(self): + Vs = FunctionSpace(self.source_mesh, "CG", 1) + Vt = FunctionSpace(self.target_mesh, "CG", 1).dual() + with self.assertRaises(ValueError) as cm: + transfer(Function(Vs), Function(Vt)) + self.assertEqual(str(cm.exception), "Spaces must be both primal or both dual.") + @parameterized.expand(["interpolate", "project"]) def test_no_sub_source_space(self, transfer_method): Vs = FunctionSpace(self.source_mesh, "CG", 1) @@ -203,7 +199,7 @@ def test_no_sub_source_space_adjoint(self, transfer_method): Vt = Vt * Vt with self.assertRaises(ValueError) as cm: transfer(Cofunction(Vs.dual()), Cofunction(Vt.dual()), transfer_method) - msg = "Source space has multiple components but target space does not." + msg = "Target space has multiple components but source space does not." self.assertEqual(str(cm.exception), msg) @parameterized.expand(["interpolate", "project"]) @@ -223,7 +219,7 @@ def test_no_sub_target_space_adjoint(self, transfer_method): Vs = Vs * Vs with self.assertRaises(ValueError) as cm: transfer(Cofunction(Vs.dual()), Cofunction(Vt.dual()), transfer_method) - msg = "Target space has multiple components but source space does not." + msg = "Source space has multiple components but target space does not." self.assertEqual(str(cm.exception), msg) @parameterized.expand(["interpolate", "project"]) @@ -255,16 +251,14 @@ def test_transfer_same_space(self, transfer_method): expected = source self.assertAlmostEqual(errornorm(expected, target), 0) - @parameterized.expand(["project"]) # TODO: interpolate (#113) + @parameterized.expand(["interpolate", "project"]) def test_transfer_same_space_adjoint(self, transfer_method): - pytest.skip() # TODO: (#114) Vs = FunctionSpace(self.source_mesh, "CG", 1) - source = Function(Vs).interpolate(self.sinusoid()) - source = function2cofunction(source) + expected = Function(Vs).interpolate(self.sinusoid()) + source = function2cofunction(expected) target = Cofunction(Vs.dual()) transfer(source, target, transfer_method) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) + self.assertAlmostEqual(errornorm(expected, cofunction2function(target)), 0) @parameterized.expand(["project", "interpolate"]) def test_transfer_same_space_mixed(self, transfer_method): @@ -279,20 +273,18 @@ def test_transfer_same_space_mixed(self, transfer_method): expected = source self.assertAlmostEqual(errornorm(expected, target), 0) - @parameterized.expand(["project"]) # TODO: interpolate (#113) + @parameterized.expand(["interpolate", "project"]) def test_transfer_same_space_mixed_adjoint(self, transfer_method): - pytest.skip() # TODO: (#114) P1 = FunctionSpace(self.source_mesh, "CG", 1) Vs = P1 * P1 - source = Function(Vs) - s1, s2 = source.subfunctions - s1.interpolate(self.sinusoid()) - s2.interpolate(-self.sinusoid()) - source = function2cofunction(source) + expected = Function(Vs) + e1, e2 = expected.subfunctions + e1.interpolate(self.sinusoid()) + e2.interpolate(-self.sinusoid()) + source = function2cofunction(expected) target = Cofunction(Vs.dual()) transfer(source, target, transfer_method) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) + self.assertAlmostEqual(errornorm(expected, cofunction2function(target)), 0) @parameterized.expand(["interpolate", "project"]) def test_transfer_same_mesh(self, transfer_method): @@ -307,7 +299,7 @@ def test_transfer_same_mesh(self, transfer_method): expected = Function(Vt).project(source) self.assertAlmostEqual(errornorm(expected, target), 0) - @parameterized.expand(["project"]) # TODO: interpolate (#113) + @parameterized.expand(["interpolate", "project"]) def test_transfer_same_mesh_adjoint(self, transfer_method): pytest.skip() # TODO: (#114) Vs = FunctionSpace(self.source_mesh, "CG", 1) @@ -343,7 +335,7 @@ def test_transfer_same_mesh_mixed(self, transfer_method): e2.project(s2) self.assertAlmostEqual(errornorm(expected, target), 0) - @parameterized.expand(["project"]) # TODO: interpolate (#113) + @parameterized.expand(["interpolate", "project"]) def test_transfer_same_mesh_mixed_adjoint(self, transfer_method): pytest.skip() # TODO: (#114) P1 = FunctionSpace(self.source_mesh, "CG", 1) @@ -377,11 +369,11 @@ def check_no_new_extrema(source, target, tol=1.0e-08): @parameterized.expand([(True, True), (True, False), (False, True), (False, False)]) def test_supermesh_project(self, same_mesh, same_degree, tol=1.0e-07): Vs = FunctionSpace(self.source_mesh, "CG", 1) + x, y = ufl.SpatialCoordinate(self.source_mesh) + source = Function(Vs).interpolate(self.sinusoid()) target_mesh = self.source_mesh if same_mesh else self.target_mesh target_degree = 1 if same_degree else 0 Vt = FunctionSpace(target_mesh, "DG", target_degree) - x, y = ufl.SpatialCoordinate(self.source_mesh) - source = Function(Vs).interpolate(self.sinusoid()) target = Function(Vt) _supermesh_project(source, target, bounded=False) expected = Function(Vt).project(source) @@ -393,10 +385,11 @@ def test_supermesh_project(self, same_mesh, same_degree, tol=1.0e-07): @parameterized.expand([(True,), (False,)]) def test_mass_lumping(self, same_mesh, tol=1.0e-08): Vs = FunctionSpace(self.source_mesh, "CG", 1) - target_mesh = self.source_mesh if same_mesh else self.target_mesh - Vt = FunctionSpace(target_mesh, "CG", 1) x, y = ufl.SpatialCoordinate(self.source_mesh) source = Function(Vs).interpolate(self.sinusoid()) - target = project(source, Vt, bounded=True) + target_mesh = self.source_mesh if same_mesh else self.target_mesh + Vt = FunctionSpace(target_mesh, "CG", 1) + target = Function(Vt) + project(source, target, bounded=True) self.assertTrue(self.check_conservation(source, target, tol=tol)) self.assertTrue(self.check_no_new_extrema(source, target, tol=tol))