From a25722249b3d2b1cdf1eea575e0b39d9144c907a Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Mon, 16 Feb 2026 15:35:11 +0000 Subject: [PATCH 01/17] Fix for Eigen compilation errors (#223) This seems to fix compilation errors that currently break the tests. --- animate/cxx/metric2d.cxx | 2 +- animate/cxx/metric3d.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/animate/cxx/metric2d.cxx b/animate/cxx/metric2d.cxx index 86d1b3e1..3e51c0aa 100644 --- a/animate/cxx/metric2d.cxx +++ b/animate/cxx/metric2d.cxx @@ -146,7 +146,7 @@ void intersect(double M_[4], const double * A_, const double * B_) { // Solve eigenvalue problem for triple product of inverse square root metric and the second metric SelfAdjointEigenSolver> eigensolver2(Sqi.transpose() * B * Sqi); Q = eigensolver2.eigenvectors(); - D = eigensolver2.eigenvalues().array().max(1).matrix().asDiagonal(); + D = eigensolver2.eigenvalues().array().cwiseMax(1).matrix().asDiagonal(); // Compute metric intersection M = Sq.transpose() * Q * D * Q.transpose() * Sq; diff --git a/animate/cxx/metric3d.cxx b/animate/cxx/metric3d.cxx index 78e0cc71..8bc60eb2 100644 --- a/animate/cxx/metric3d.cxx +++ b/animate/cxx/metric3d.cxx @@ -187,7 +187,7 @@ void intersect(double M_[9], const double * A_, const double * B_) { // Solve eigenvalue problem for triple product of inverse square root metric and the second metric SelfAdjointEigenSolver> eigensolver2(Sqi.transpose() * B * Sqi); Q = eigensolver2.eigenvectors(); - D = eigensolver2.eigenvalues().array().max(1).matrix().asDiagonal(); + D = eigensolver2.eigenvalues().array().cwiseMax(1).matrix().asDiagonal(); // Compute metric intersection M = Sq.transpose() * Q * D * Q.transpose() * Sq; From da286015480362ebddd5e9da5fac8120ddc8bd72 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+joewallwork@users.noreply.github.com> Date: Mon, 23 Feb 2026 11:41:58 +0000 Subject: [PATCH 02/17] Separate `main` and `release` workflows (#218) See https://github.com/mesh-adaptation/docs/issues/156. This PR adds a separate test suite workflow for the new `release` branch. That test suite not only runs on Animate's `release` branch, but also with the `release` version of Firedrake. --- .../{test_suite.yml => test_suite_main.yml} | 8 +++---- .github/workflows/test_suite_release.yml | 24 +++++++++++++++++++ 2 files changed, 28 insertions(+), 4 deletions(-) rename .github/workflows/{test_suite.yml => test_suite_main.yml} (71%) create mode 100644 .github/workflows/test_suite_release.yml diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite_main.yml similarity index 71% rename from .github/workflows/test_suite.yml rename to .github/workflows/test_suite_main.yml index d4f4ef5a..83299516 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite_main.yml @@ -1,14 +1,14 @@ -name: 'Run Animate Test Suite' +name: 'Run Animate test suite (main)' on: - # Run test suite whenever main is updated + # Run test suite whenever the main branch is updated push: branches: - main paths: - '**.py' - '**.cxx' - - '.github/workflows/test_suite.yml' + - '.github/workflows/test_suite_main.yml' - '.gitmodules' - 'pyproject.toml' @@ -17,7 +17,7 @@ on: paths: - '**.py' - '**.cxx' - - '.github/workflows/test_suite.yml' + - '.github/workflows/test_suite_main.yml' - '.gitmodules' - 'pyproject.toml' diff --git a/.github/workflows/test_suite_release.yml b/.github/workflows/test_suite_release.yml new file mode 100644 index 00000000..af664a98 --- /dev/null +++ b/.github/workflows/test_suite_release.yml @@ -0,0 +1,24 @@ +name: 'Run Animate test suite (release)' + +on: + # Run test suite whenever the release branch is updated + push: + branches: + - release + paths: + - '.github/workflows/test_suite_release.yml' + + # Run test suite whenever commits that change this workflow are pushed to an open PR + pull_request: + paths: + - '.github/workflows/test_suite_release.yml' + + # Run test suite at 3AM on the first of the month + schedule: + - cron: '0 3 1 * *' + +jobs: + test_suite: + uses: mesh-adaptation/docs/.github/workflows/reusable_test_suite.yml@main + with: + docker-image: 'firedrake-parmmg:release' From 30010ee2a4b048898f78047c7f8a5f296f535d8a Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Mon, 23 Feb 2026 17:18:36 +0000 Subject: [PATCH 03/17] Fix import of cached_property (#229) Due to https://github.com/firedrakeproject/firedrake/pull/4903 --- animate/adapt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/animate/adapt.py b/animate/adapt.py index 7bd25c82..5a8bde63 100644 --- a/animate/adapt.py +++ b/animate/adapt.py @@ -1,11 +1,11 @@ import abc import os +from functools import cached_property from shutil import rmtree import firedrake.checkpointing as fchk import firedrake.functionspace as ffs import firedrake.mesh as fmesh -import firedrake.utils as futils from firedrake import COMM_SELF, COMM_WORLD from firedrake.cython.dmcommon import to_petsc_local_numbering from firedrake.petsc import PETSc @@ -87,7 +87,7 @@ def __init__(self, mesh, metric, name=None, comm=None): self.metric = metric self.projectors = [] - @futils.cached_property + @cached_property @PETSc.Log.EventDecorator() def adapted_mesh(self): """ From 98ec8109de7b474162063a3deaf8e88b8de4ab67 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+joewallwork@users.noreply.github.com> Date: Tue, 24 Feb 2026 15:28:43 +0000 Subject: [PATCH 04/17] Move `to_petsc_local_numbering` into Animate (#226) Closes #225. Moves `to_petsc_local_numbering` to Animate and adds a basic unit test. --- .gitignore | 3 ++ animate/adapt.py | 2 +- animate/cython/numbering.pyx | 59 +++++++++++++++++++++++++ pyproject.toml | 34 ++++++++++++--- setup.py | 83 ++++++++++++++++++++++++++++++++++++ test/test_numbering.py | 47 ++++++++++++++++++++ 6 files changed, 222 insertions(+), 6 deletions(-) create mode 100644 animate/cython/numbering.pyx create mode 100644 setup.py create mode 100644 test/test_numbering.py diff --git a/.gitignore b/.gitignore index 66b6cac1..3e77c7bd 100644 --- a/.gitignore +++ b/.gitignore @@ -10,9 +10,12 @@ __pycache__/ .coverage # File extensions +*.c +*.h5 *.jpg *.jpeg *.pdf *.png *.pvd +*.so *.vtu diff --git a/animate/adapt.py b/animate/adapt.py index 5a8bde63..b3baa595 100644 --- a/animate/adapt.py +++ b/animate/adapt.py @@ -7,11 +7,11 @@ import firedrake.functionspace as ffs import firedrake.mesh as fmesh from firedrake import COMM_SELF, COMM_WORLD -from firedrake.cython.dmcommon import to_petsc_local_numbering from firedrake.petsc import PETSc from firedrake.projection import Projector from .checkpointing import get_checkpoint_dir, load_checkpoint, save_checkpoint +from .cython.numbering import to_petsc_local_numbering from .metric import RiemannianMetric __all__ = ["MetricBasedAdaptor", "adapt"] diff --git a/animate/cython/numbering.pyx b/animate/cython/numbering.pyx new file mode 100644 index 00000000..3abaee2b --- /dev/null +++ b/animate/cython/numbering.pyx @@ -0,0 +1,59 @@ +""" +Cython module for handling numberings of PETSc Vecs corresponding to Firedrake Functions. +""" +import numpy as np +from firedrake.petsc import PETSc +from firedrake.utils import IntType, ScalarType + +cimport numpy as np +cimport petsc4py.PETSc as PETSc + +from petsc4py.PETSc cimport CHKERR + + +cdef extern from "petsc.h": + ctypedef long PetscInt + ctypedef enum PetscErrorCode: + PETSC_SUCCESS + PETSC_ERR_LIB + + +cdef extern from "petscis.h" nogil: + PetscErrorCode PetscSectionGetDof(PETSc.PetscSection,PetscInt,PetscInt*) + PetscErrorCode PetscSectionGetOffset(PETSc.PetscSection,PetscInt,PetscInt*) + + +def to_petsc_local_numbering(PETSc.Vec vec, V): + """ + Reorder a PETSc Vec corresponding to a Firedrake Function w.r.t. + the PETSc natural numbering. + + :arg vec: the PETSc Vec to reorder; must be a global vector + :arg V: the FunctionSpace of the Function which the Vec comes from + :ret out: a copy of the Vec, ordered with the PETSc natural numbering + """ + cdef int dim, idx, start, end, p, d, k + cdef PetscInt dof, off + cdef PETSc.Vec out + cdef PETSc.Section section + cdef np.ndarray varray, oarray + + section = V.dm.getGlobalSection() + out = vec.duplicate() + varray = vec.array_r + oarray = out.array + dim = V.value_size + idx = 0 + start, end = vec.getOwnershipRange() + for p in range(*section.getChart()): + CHKERR(PetscSectionGetDof(section.sec, p, &dof)) + if dof > 0: + CHKERR(PetscSectionGetOffset(section.sec, p, &off)) + assert off >= 0 + off *= dim + for d in range(dof): + for k in range(dim): + oarray[idx] = varray[off + dim * d + k - start] + idx += 1 + assert idx == (end - start) + return out diff --git a/pyproject.toml b/pyproject.toml index 581e4e99..8f3e0d11 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,18 @@ [build-system] -requires = ["setuptools>=77.0.3"] +build-backend = "setuptools.build_meta" +requires = [ + "Cython", + "firedrake", + "numpy", + "petsc4py", + "petsctools", + "setuptools>=77.0.3", + "wheel", +] [project] name = "animate" -version = "0.3" +version = "0.4" authors = [ {name = "Joseph G. Wallwork", email = "joe.wallwork@outlook.com"}, {name = "Davor Dundovic"}, @@ -43,11 +52,25 @@ Homepage = "https://mesh-adaptation.github.io" Documentation = "https://mesh-adaptation.github.io/animate/index.html" Repository = "https://github.com/mesh-adaptation/animate" +# Note: +# - Extension modules are generated dynamically at build time using the legacy +# helper (`setup.py`) which performs PETSc detection and populates Extension +# objects based on the local PETSc installation (via petsctools/petsc4py). +# - We intentionally do not declare ext-modules statically here so the build +# helper can compute include/library dirs and other flags based on an +# existing PETSc build. +# +# The build backend requirements above ensure that, during isolated builds, +# the build environment will have Cython, numpy, petsctools, petsc4py, etc. +# available before the build helper runs Cython to generate .c files. [tool.setuptools] packages = ["animate", "adapt_common"] [tool.setuptools.package-data] -animate = ["cxx/*.cxx"] +animate = [ + "cxx/*.cxx", + "cython/*.pyx", +] [tool.setuptools.package-dir] adapt_common = "adapt_common/adapt_common" @@ -63,6 +86,7 @@ select = [ "F", # Pyflakes "I", # isort ] + [tool.ruff.lint.per-file-ignores] "demos/*" = [ "E402", # module level import not at top of file @@ -72,6 +96,6 @@ select = [ [tool.pytest.ini_options] filterwarnings = [ - "ignore:`np.bool8` is a deprecated alias for `np.bool_`*:DeprecationWarning", - "ignore:unable to find git revision*:UserWarning", + "ignore:`np.bool8` is a deprecated alias for `np.bool_`*:DeprecationWarning", + "ignore:unable to find git revision*:UserWarning", ] diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..bfd6d59f --- /dev/null +++ b/setup.py @@ -0,0 +1,83 @@ +"""Setup script for building Cython extensions for the `animate` package. + +The script performs PETSc detection and populates Extension objects based on the local +PETSc installation (via petsctools/petsc4py). +""" + +import os +from dataclasses import dataclass, field + +import numpy as np +import petsc4py +import petsctools +from Cython.Build import cythonize +from setuptools import Extension, find_packages, setup + + +@dataclass +class ExternalDependency: + """Courtesy of Firedrake. See + https://github.com/firedrakeproject/firedrake/blob/main/setup.py + """ + + include_dirs: list[str] = field(default_factory=list, init=True) + extra_compile_args: list[str] = field(default_factory=list, init=True) + libraries: list[str] = field(default_factory=list, init=True) + library_dirs: list[str] = field(default_factory=list, init=True) + extra_link_args: list[str] = field(default_factory=list, init=True) + runtime_library_dirs: list[str] = field(default_factory=list, init=True) + + def __add__(self, other): + combined = {} + for f in self.__dataclass_fields__.keys(): + combined[f] = getattr(self, f) + getattr(other, f) + return self.__class__(**combined) + + def keys(self): + return self.__dataclass_fields__.keys() + + def __getitem__(self, key): + try: + return getattr(self, key) + except AttributeError as attr_err: + raise KeyError(f"Key {key} not present") from attr_err + + +def extensions(): + """Returns a list of Cython extensions to be compiled.""" + + # Define external dependencies + mpi_ = ExternalDependency( + extra_compile_args=petsctools.get_petscvariables()["MPICC_SHOW"].split()[1:], + ) + numpy_ = ExternalDependency(include_dirs=[np.get_include()]) + petsc_dir = petsctools.get_petsc_dir() + petsc_dirs = [petsc_dir, os.path.join(petsc_dir, petsctools.get_petsc_arch())] + petsc_includes = [petsc4py.get_include()] + [ + os.path.join(d, "include") for d in petsc_dirs + ] + petsc_ = ExternalDependency( + libraries=["petsc"], + include_dirs=petsc_includes, + library_dirs=[os.path.join(petsc_dirs[-1], "lib")], + runtime_library_dirs=[os.path.join(petsc_dirs[-1], "lib")], + ) + + # Define Cython extensions + cython_list = [ + Extension( + name="animate.cython.numbering", + language="c", + sources=[os.path.join("animate", "cython", "numbering.pyx")], + **(mpi_ + petsc_ + numpy_), + ) + ] + return cythonize(cython_list) + + +if __name__ == "__main__": + # Run the setup function to build the Cython extensions + setup( + packages=find_packages(), + ext_modules=extensions(), + ) diff --git a/test/test_numbering.py b/test/test_numbering.py new file mode 100644 index 00000000..9d859c30 --- /dev/null +++ b/test/test_numbering.py @@ -0,0 +1,47 @@ +"""Module containing unit tests for the cython.numbering module.""" + +import firedrake as fd +import numpy as np +import pytest +from animate.cython.numbering import to_petsc_local_numbering + + +@pytest.fixture +def mesh(): + """Fixture to create a uniform unit square mesh.""" + return fd.UnitSquareMesh(4, 4) + + +@pytest.fixture(params=[0, 1, 2]) +def rank(request): + """Fixture specifying function space rank.""" + return request.param + + +@pytest.fixture +def function_space(mesh, rank): + """Fixture to create a function space on a given mesh with a given rank.""" + try: + return { + 0: fd.FunctionSpace, + 1: fd.VectorFunctionSpace, + 2: fd.TensorFunctionSpace, + }[rank](mesh, "CG", 1) + except KeyError as key_err: + raise ValueError(f"Rank {rank} not considered.") from key_err + + +def test_is_permutation(function_space): + """Test that to_petsc_local_numbering is indeed a permutation.""" + f = fd.Function(function_space) + + # Fill the vector with arbitrary values + with f.dat.vec as vec: + vec.array[:] = np.arange(vec.size) + + # Verify that reordering according to the PETSc numbering is a permutation + with f.dat.vec_ro as vec: + reordered_vec = to_petsc_local_numbering(vec, function_space) + assert reordered_vec is not None + assert reordered_vec.size == vec.size + assert sorted(reordered_vec.array) == sorted(vec.array) From a9300ff27a098c5597a631c3a35cbd5d4085a549 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Mon, 2 Mar 2026 15:11:28 +0000 Subject: [PATCH 05/17] Provide parmmg with complete local vector including halos (#228) This changes to_petsc_local_numbering to take in a local vector (i.e. a sequential vector that includes halo DOFs) and returns a local vector in the "natural ordering" of the dmplex topological points. This is what the dmplex parmmg interface expects from the metric input vector. Should not change anything in serial. Appears to fix the warnings about indefinite metric in parmmg tests - which I think were caused by halo entries that were copied, out-of-bounds, from the global vector we were providing previously Closes #181 Partial addresses #215 The change to to_petsc_local_numbering and the added parallel tests, seems to make the stalls in the CI (#215) worse: happening most CI runs, rather than 50% of the time. Therefore changes in the (unrelated) serial checkpointing code have been included, which _seem_ to make the CI pass more reliably (not failing so far on ~8 runs) - but it is quite possible (likely?) that the underlying issue is not resolved. --- animate/adapt.py | 13 +++++++--- animate/cython/numbering.pyx | 21 ++++++++++------ test/test_numbering.py | 47 +++++++++++++++++++++++++++++++----- 3 files changed, 65 insertions(+), 16 deletions(-) diff --git a/animate/adapt.py b/animate/adapt.py index b3baa595..e9ed3007 100644 --- a/animate/adapt.py +++ b/animate/adapt.py @@ -1,4 +1,5 @@ import abc +import gc import os from functools import cached_property from shutil import rmtree @@ -97,10 +98,9 @@ def adapted_mesh(self): :rtype: :class:`firedrake.mesh.MeshGeometry` """ self.metric.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) - size = self.metric.dat.dataset.layout_vec.getSizes() - data = self.metric.dat._data[: size[0]] + data = self.metric.dat.data_ro_with_halos v = PETSc.Vec().createWithArray( - data, size=size, bsize=self.metric.dat.cdim, comm=self.mesh.comm + data, size=data.size, bsize=self.metric.dat.cdim, comm=COMM_SELF ) reordered = to_petsc_local_numbering(v, self.metric.function_space()) v.destroy() @@ -193,13 +193,20 @@ def adapt(mesh, *metrics, name=None, serialise=None, remove_checkpoints=True): chk_fpath = os.path.join(chk_dir, "adapted_mesh_checkpoint.h5") metric_name = "tmp_metric" save_checkpoint(chk_fpath, metric, metric_name) + # Ensure all processes are finished writing + COMM_WORLD.barrier() if COMM_WORLD.rank == 0: metric0 = load_checkpoint(chk_fpath, mesh.name, metric_name, comm=COMM_SELF) adaptor0 = MetricBasedAdaptor(metric0._mesh, metric0, name=name) with fchk.CheckpointFile(chk_fpath, "w", comm=COMM_SELF) as chk: chk.save_mesh(adaptor0.adapted_mesh) + # Ensure rank 0 is finished writing COMM_WORLD.barrier() + # Garbage collection might be called at different times on different ranks due + # to diverging paths, which appears to stall in final cleanup on Python + # system exit. Ensure everything is in-sync again at this point + gc.collect() # In parallel, load from the checkpoint if not os.path.exists(chk_fpath): diff --git a/animate/cython/numbering.pyx b/animate/cython/numbering.pyx index 3abaee2b..c80b1248 100644 --- a/animate/cython/numbering.pyx +++ b/animate/cython/numbering.pyx @@ -26,25 +26,27 @@ cdef extern from "petscis.h" nogil: def to_petsc_local_numbering(PETSc.Vec vec, V): """ Reorder a PETSc Vec corresponding to a Firedrake Function w.r.t. - the PETSc natural numbering. + the PETSc natural numbering, i.e. the numbering consistent with + that of the DMPlex topological points. - :arg vec: the PETSc Vec to reorder; must be a global vector + :arg vec: the PETSc Vec to reorder; must be a local vector, i.e. + a sequential vector that includes all (owned and halo) DoFs :arg V: the FunctionSpace of the Function which the Vec comes from :ret out: a copy of the Vec, ordered with the PETSc natural numbering """ - cdef int dim, idx, start, end, p, d, k + cdef int dim, idx, lsize, p, d, k cdef PetscInt dof, off cdef PETSc.Vec out cdef PETSc.Section section cdef np.ndarray varray, oarray - section = V.dm.getGlobalSection() + section = V.dm.getLocalSection() out = vec.duplicate() varray = vec.array_r oarray = out.array dim = V.value_size idx = 0 - start, end = vec.getOwnershipRange() + lsize = vec.getSize() for p in range(*section.getChart()): CHKERR(PetscSectionGetDof(section.sec, p, &dof)) if dof > 0: @@ -53,7 +55,12 @@ def to_petsc_local_numbering(PETSc.Vec vec, V): off *= dim for d in range(dof): for k in range(dim): - oarray[idx] = varray[off + dim * d + k - start] + oarray[idx] = varray[off + dim * d + k] idx += 1 - assert idx == (end - start) + if idx != lsize: + raise ValueError( + f"Number of local section entries not the same as vector size" + f"({idx} vs. {lsize}). Need to provide local vector including halo DoFs." + ) + return out diff --git a/test/test_numbering.py b/test/test_numbering.py index 9d859c30..130411a2 100644 --- a/test/test_numbering.py +++ b/test/test_numbering.py @@ -4,6 +4,7 @@ import numpy as np import pytest from animate.cython.numbering import to_petsc_local_numbering +from firedrake.petsc import PETSc @pytest.fixture @@ -35,13 +36,47 @@ def test_is_permutation(function_space): """Test that to_petsc_local_numbering is indeed a permutation.""" f = fd.Function(function_space) + comm = function_space.comm + # unique float 0<=x<1 for each rank: + rank_fraction = comm.rank / (comm.size + 1) + + # Initialise owned and in particular halo DoFs, so we can check they + # have been updated later on + f.assign(-1) + # Fill the vector with arbitrary values + # f.dat.vec below is a global Vec, so owned entries only + # halos values should be updated automatically coming out of the context with f.dat.vec as vec: - vec.array[:] = np.arange(vec.size) + vec.array[:] = np.arange(vec.sizes[0]) + rank_fraction + + # Check that this is the case: + owned = f.dat.data_ro.flatten() + halos = f.dat.data_ro_with_halos[owned.size:].flatten() + np.testing.assert_equal(owned, np.arange(owned.size) + rank_fraction) + if halos.size > 0: + # all halo enties should be set and come from different rank + assert all(halos >= 0.0) + assert all(np.modf(halos)[0] != rank_fraction) # Verify that reordering according to the PETSc numbering is a permutation - with f.dat.vec_ro as vec: - reordered_vec = to_petsc_local_numbering(vec, function_space) - assert reordered_vec is not None - assert reordered_vec.size == vec.size - assert sorted(reordered_vec.array) == sorted(vec.array) + data = f.dat.data_ro_with_halos + lvec = PETSc.Vec().createWithArray( + data, size=data.size, bsize=f.dat.cdim, comm=fd.COMM_SELF + ) + reordered_vec = to_petsc_local_numbering(lvec, function_space) + assert reordered_vec is not None + assert reordered_vec.size == lvec.size + assert sorted(reordered_vec.array) == sorted(lvec.array) + lvec.destroy() + reordered_vec.destroy() + + +@pytest.mark.parallel(nprocs=2) +def test_is_permutation_np2(function_space): + test_is_permutation(function_space) + + +@pytest.mark.parallel(nprocs=3) +def test_is_permutation_np3(function_space): + test_is_permutation(function_space) From 452438fc1b457d930e71d51be4830d9ed4456b1d Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Tue, 3 Mar 2026 08:59:52 +0000 Subject: [PATCH 06/17] Raise an error if field isn't scalar and test --- animate/metric.py | 8 ++++++++ test/test_metric.py | 28 +++++++++++++++++++++++++++- 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/animate/metric.py b/animate/metric.py index 2119fc6e..a90c4b0c 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -349,6 +349,14 @@ def compute_hessian(self, field, method="mixed_L2", **kwargs): for the gradient recovery. The target space for the Hessian recovery is inherited from the metric itself. """ + if len(field.function_space().value_shape) > 0: + value_error = ( + "RiemannianMetric.compute_hessian only accepts scalar fields. To " + "recover Hessians of higher rank fields, call the method on separate " + "RiemannianMetrics for each component and combine them with " + "RiemannianMetric.combine." + ) + raise ValueError(value_error) if method == "L2": gradient = recover_gradient_l2( field, target_space=kwargs.get("target_space") diff --git a/test/test_metric.py b/test/test_metric.py index 9faeaa34..4ae2f2a9 100644 --- a/test/test_metric.py +++ b/test/test_metric.py @@ -2,7 +2,6 @@ import numpy as np import ufl -from adapt_common.reduction import function_data_min from firedrake.bcs import DirichletBC from firedrake.constant import Constant from firedrake.function import Function @@ -16,6 +15,7 @@ from sensors import bowl, hyperbolic, interweaved, multiscale from test_setup import uniform_mesh, uniform_metric +from adapt_common.reduction import function_data_min from animate.metric import P0Metric, RiemannianMetric @@ -194,6 +194,32 @@ class TestHessianMetric(MetricTestCase): Unit tests for the :meth:`compute_hessian` method of :class:`RiemannianMetric`. """ + def test_vector_rank_error(self): + mesh = uniform_mesh(2) + P1_vec = VectorFunctionSpace(mesh, "CG", 1) + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + with self.assertRaises(ValueError) as cm: + RiemannianMetric(P1_ten).compute_hessian(Function(P1_vec)) + msg = ( + "RiemannianMetric.compute_hessian only accepts scalar fields. To " + "recover Hessians of higher rank fields, call the method on separate " + "RiemannianMetrics for each component and combine them with " + "RiemannianMetric.combine." + ) + self.assertEqual(str(cm.exception), msg) + + def test_tensor_rank_error(self): + P1_ten = TensorFunctionSpace(uniform_mesh(2), "CG", 1) + with self.assertRaises(ValueError) as cm: + RiemannianMetric(P1_ten).compute_hessian(Function(P1_ten)) + msg = ( + "RiemannianMetric.compute_hessian only accepts scalar fields. To " + "recover Hessians of higher rank fields, call the method on separate " + "RiemannianMetrics for each component and combine them with " + "RiemannianMetric.combine." + ) + self.assertEqual(str(cm.exception), msg) + def test_bowl(self, dim=2, places=7): mesh = uniform_mesh(dim, 4, recentre=True) P1_ten = TensorFunctionSpace(mesh, "CG", 1) From 7a96e796fac3abdd0e45be577c274a649e8bd06a Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Thu, 5 Mar 2026 08:21:40 +0000 Subject: [PATCH 07/17] Lint --- test/test_metric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_metric.py b/test/test_metric.py index 4ae2f2a9..6dbaccd0 100644 --- a/test/test_metric.py +++ b/test/test_metric.py @@ -2,6 +2,7 @@ import numpy as np import ufl +from adapt_common.reduction import function_data_min from firedrake.bcs import DirichletBC from firedrake.constant import Constant from firedrake.function import Function @@ -15,7 +16,6 @@ from sensors import bowl, hyperbolic, interweaved, multiscale from test_setup import uniform_mesh, uniform_metric -from adapt_common.reduction import function_data_min from animate.metric import P0Metric, RiemannianMetric From ed14970614206a3ea4fd3209cf8f032df101ac80 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+joewallwork@users.noreply.github.com> Date: Fri, 6 Mar 2026 08:05:01 +0000 Subject: [PATCH 08/17] Apply suggestion from @stephankramer code review Co-authored-by: Stephan Kramer --- animate/metric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/animate/metric.py b/animate/metric.py index a90c4b0c..65f4ca9c 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -349,7 +349,7 @@ def compute_hessian(self, field, method="mixed_L2", **kwargs): for the gradient recovery. The target space for the Hessian recovery is inherited from the metric itself. """ - if len(field.function_space().value_shape) > 0: + if len(field.ufl_shape) > 0: value_error = ( "RiemannianMetric.compute_hessian only accepts scalar fields. To " "recover Hessians of higher rank fields, call the method on separate " From dcb3e165bfbd25b5e9740d6a9f7020b2d06ab5bf Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 6 Mar 2026 08:23:11 +0000 Subject: [PATCH 09/17] Fix recovery test --- test/test_recovery.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_recovery.py b/test/test_recovery.py index 1e7ebe4e..57a7fc31 100644 --- a/test/test_recovery.py +++ b/test/test_recovery.py @@ -47,7 +47,7 @@ def test_clement_function_error(self): self.assertEqual(str(cm.exception), msg) def test_clement_space_error(self): - f = self.get_func_ones("RT", 1) + f = self.get_func_ones("DG", 1) with self.assertRaises(ValueError) as cm: self.metric.compute_hessian(f, method="Clement") msg = ( From 640b04c337b71a87a47689b13b466a570be2701e Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+joewallwork@users.noreply.github.com> Date: Fri, 6 Mar 2026 11:33:53 +0000 Subject: [PATCH 10/17] Fixes for `TestRecoverySetup` (#232) Follow up for the remaining fixes after accidentally merging #231. --- test/test_recovery.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_recovery.py b/test/test_recovery.py index 57a7fc31..5c7bf619 100644 --- a/test/test_recovery.py +++ b/test/test_recovery.py @@ -15,6 +15,7 @@ from animate.math import construct_basis from animate.metric import RiemannianMetric +from animate.recovery import recover_gradient_l2 # --------------------------- # standard tests for pytest @@ -47,7 +48,7 @@ def test_clement_function_error(self): self.assertEqual(str(cm.exception), msg) def test_clement_space_error(self): - f = self.get_func_ones("DG", 1) + f = self.get_func_ones("DG", 0) with self.assertRaises(ValueError) as cm: self.metric.compute_hessian(f, method="Clement") msg = ( @@ -69,14 +70,14 @@ def test_clement_degree_error(self): def test_l2_projection_function_error(self): f = bowl(*ufl.SpatialCoordinate(self.mesh)) with self.assertRaises(ValueError) as cm: - self.metric.compute_hessian(f, method="L2") + recover_gradient_l2(f) msg = "If a target space is not provided then the input must be a Function." self.assertEqual(str(cm.exception), msg) def test_l2_projection_rank_error(self): f = Function(TensorFunctionSpace(self.mesh, "CG", 1)) with self.assertRaises(ValueError) as cm: - self.metric.compute_hessian(f, method="L2") + recover_gradient_l2(f) msg = ( "L2 projection can only be used to compute gradients of scalar or vector" " Functions, not Functions of rank 2." From eadf7fd8ba27175a4c071d40209ff26b6669aa4e Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Fri, 6 Mar 2026 14:42:11 +0000 Subject: [PATCH 11/17] Fix intersection of multiple metrics (#235) Closes #234 Make the existing test fail on this bug by reordering, and add a test with multiple spatially varying metrics. --- animate/metric.py | 37 ++++++++++++++++++++----------------- test/test_metric.py | 27 ++++++++++++++++++++++++--- 2 files changed, 44 insertions(+), 20 deletions(-) diff --git a/animate/metric.py b/animate/metric.py index 65f4ca9c..b449dbd1 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -629,25 +629,28 @@ def intersect(self, *metrics): "Cannot intersect metrics with different function spaces." ) - # Intersect the metrics recursively one at a time if len(metrics) == 0: - pass - elif len(metrics) == 1: - v1 = self._create_from_array(self.dat.data_with_halos) - v2 = self._create_from_array(metrics[0].dat.data_ro_with_halos) - vout = self._create_from_array(np.zeros_like(self.dat.data_with_halos)) - - # Compute the intersection on the PETSc level - self._plex.metricIntersection2(v1, v2, vout) - - # Assign to the output of the intersection - size = np.shape(self.dat.data_with_halos) - self.dat.data_with_halos[:] = np.reshape(vout.array, size) - v2.destroy() - v1.destroy() - vout.destroy() - else: + return self + + # Intersect the metrics recursively, starting with metrics[0] + v1 = self._create_from_array(self.dat.data_with_halos) + v2 = self._create_from_array(metrics[0].dat.data_ro_with_halos) + vout = self._create_from_array(np.zeros_like(self.dat.data_with_halos)) + + # Compute the intersection on the PETSc level + self._plex.metricIntersection2(v1, v2, vout) + + # Assign to the output of the intersection + size = np.shape(self.dat.data_with_halos) + self.dat.data_with_halos[:] = np.reshape(vout.array, size) + v2.destroy() + v1.destroy() + vout.destroy() + + # Intersect with the remaining metrics + if len(metrics) > 1: self.intersect(*metrics[1:]) + return self @PETSc.Log.EventDecorator() diff --git a/test/test_metric.py b/test/test_metric.py index 6dbaccd0..063f60e3 100644 --- a/test/test_metric.py +++ b/test/test_metric.py @@ -308,16 +308,37 @@ def test_multiple_intersect(self, dim): mesh = uniform_mesh(dim, 1) P1_ten = TensorFunctionSpace(mesh, "CG", 1) - metric1 = uniform_metric(P1_ten, 100.0) - metric2 = uniform_metric(P1_ten, 40.0) + metric1 = uniform_metric(P1_ten, 40.0) + metric2 = uniform_metric(P1_ten, 100.0) metric3 = uniform_metric(P1_ten, 20.0) - expected = metric1 + expected = metric2 metric = RiemannianMetric(P1_ten) metric.assign(metric1) metric.intersect(metric2, metric3) self.assertAlmostMatching(metric, expected) + @parameterized.expand([[2], [3]]) + def test_multiple_variable_intersect(self, dim): + mesh = uniform_mesh(dim, 1) + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + + expected = RiemannianMetric(P1_ten) + expected.interpolate(ufl.Identity(dim)) + metrics = [] + for i in range(P1_ten.node_count): + metric = RiemannianMetric(P1_ten) + metric.interpolate(ufl.Identity(dim)) + # change to (i+1)I at node i only + metric.dat.data[i] *= float(i+1) + expected.dat.data[i] *= float(i+1) + metrics.append(metric) + + metric = RiemannianMetric(P1_ten) + metric.assign(metrics[0]) + metric.intersect(*metrics[1:]) + self.assertAlmostMatching(metric, expected) + class TestNormalisation(MetricTestCase): """ From 59e4da13607d1660766d84e06ec6542c44e5c6d1 Mon Sep 17 00:00:00 2001 From: Sia Ghelichkhan Date: Thu, 2 Apr 2026 22:37:32 +1100 Subject: [PATCH 12/17] Use petsctools.get_petsc_dirs instead of pyop2.utils.get_petsc_dir Firedrake PR firedrakeproject/firedrake#5002 removed get_petsc_dir from pyop2.utils in favour of the standalone petsctools package. The setup.py in this repo already uses petsctools, but quality.py was still importing from pyop2.utils, which now breaks at import time. Fixes g-adopt/g-adopt#492 Co-authored-by: Stephan Kramer --- animate/quality.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/animate/quality.py b/animate/quality.py index d3e060f4..2cd6e298 100644 --- a/animate/quality.py +++ b/animate/quality.py @@ -8,10 +8,10 @@ import ufl from firedrake.__future__ import interpolate from firedrake.petsc import PETSc +from petsctools import get_petsc_dirs from pyop2 import op2 -from pyop2.utils import get_petsc_dir -petsc_dirs = get_petsc_dir() +petsc_dirs = get_petsc_dirs() include_dir = ["%s/include/eigen3" % petsc_dirs[-1]] __all__ = ["QualityMeasure"] From ac3aeff3b455c52e4460b3030eaa96d0dbe560f7 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Wed, 8 Apr 2026 16:47:10 +0100 Subject: [PATCH 13/17] Change to reflect new checkpoint layout Apparently topologies are now stored under their topology dm name. No idea why we're actually checking this.... --- test/test_checkpoint.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_checkpoint.py b/test/test_checkpoint.py index f687f356..9f4e1109 100644 --- a/test/test_checkpoint.py +++ b/test/test_checkpoint.py @@ -7,7 +7,6 @@ import ufl from firedrake.function import Function from firedrake.functionspace import FunctionSpace, TensorFunctionSpace -from firedrake.mesh import _generate_default_mesh_topology_name from test_setup import uniform_mesh from animate.checkpointing import get_checkpoint_dir, load_checkpoint, save_checkpoint @@ -57,8 +56,8 @@ def test_save(self, filename="test_save.h5", metric=None): metric = metric or self.metric save_checkpoint(fpath, metric) self.assertTrue(os.path.exists(fpath)) - mesh_name = metric._mesh.name - topology_name = _generate_default_mesh_topology_name(mesh_name) + mesh = ufl.domain.extract_unique_domain(metric) + topology_name = mesh.topology_dm.name with h5py.File(fpath, "r") as h5: self.assertTrue("topologies" in h5) self.assertTrue(topology_name in h5["topologies"].keys()) From 035de7f75189446f924dbc10445ded46e734c6ac Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Wed, 8 Apr 2026 17:21:25 +0100 Subject: [PATCH 14/17] Loosen tolerance Precision here hinges on the default solver parameters used in _compute_gradient_and_hessian which currently falls back on firedrake default for tolerance in outer iteration (rtol=1e-7). This appears to give an error ~ 8e-8, which is rounded to 1e-7 and thus fails the assert almost equal test with places=7. Setting a tighter tolerance indeed reduces the error as expected. Not entirely sure why the error has increased, but presumably was close previously as well but <5e-8 and thus passing the test. --- test/test_metric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_metric.py b/test/test_metric.py index 063f60e3..e0b3bf72 100644 --- a/test/test_metric.py +++ b/test/test_metric.py @@ -220,7 +220,7 @@ def test_tensor_rank_error(self): ) self.assertEqual(str(cm.exception), msg) - def test_bowl(self, dim=2, places=7): + def test_bowl(self, dim=2, places=6): mesh = uniform_mesh(dim, 4, recentre=True) P1_ten = TensorFunctionSpace(mesh, "CG", 1) metric = RiemannianMetric(P1_ten).compute_hessian(bowl(*mesh.coordinates)) From 72d13d517413da96165caa5311ec0e2a91248cc4 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Thu, 11 Jun 2026 15:37:04 +0100 Subject: [PATCH 15/17] Introduce separate requirements-build.txt (#241) Addresses #240 This is needed because we are instructing our users to build with `pip install --no-build-isolation` which tells pip to build in the current python env and ignores the build dependencies specified in pyproject.toml. Following the installation instructions for Firedrake release however, this environment may not have these build requirements installed already (notably Cython and a sufficiently new setuptools). Therefore we need to instruct users to run `pip install -r animate/requirements-build.txt` beforehand. TODO (after merge): * update wiki installation instructions * use the additional step in the CI --- requirements-build.txt | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 requirements-build.txt diff --git a/requirements-build.txt b/requirements-build.txt new file mode 100644 index 00000000..4cae0098 --- /dev/null +++ b/requirements-build.txt @@ -0,0 +1,11 @@ +# Build requirements needed to build cython module. +# Note that we assume we are in a valid Firedrake +# environment that have firedrake and petsc4py installed +# already - we do not specify these here as they may +# cause the locally installed ones to be replaced by the +# latest release from PyPI +Cython +numpy +petsctools +setuptools>=77.0.3 +wheel From 531d51ca371c398d1ec98f8e0c5a90a7f4195d4c Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Sat, 13 Jun 2026 14:37:16 +0100 Subject: [PATCH 16/17] Revert "Introduce separate requirements-build.txt (#242)" This reverts commit 6c141066c1c2610b979e6e61d9749c150c4be788. --- requirements-build.txt | 11 ----------- 1 file changed, 11 deletions(-) delete mode 100644 requirements-build.txt diff --git a/requirements-build.txt b/requirements-build.txt deleted file mode 100644 index 4cae0098..00000000 --- a/requirements-build.txt +++ /dev/null @@ -1,11 +0,0 @@ -# Build requirements needed to build cython module. -# Note that we assume we are in a valid Firedrake -# environment that have firedrake and petsc4py installed -# already - we do not specify these here as they may -# cause the locally installed ones to be replaced by the -# latest release from PyPI -Cython -numpy -petsctools -setuptools>=77.0.3 -wheel From 20513d9fe1074044e35a6e9fa37d7ed5382a07e8 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Sat, 13 Jun 2026 14:37:20 +0100 Subject: [PATCH 17/17] Revert "Fix release branch (#230)" This reverts commit 68242c8855ed74b9198b85c28a6d110ce03645a4. --- .../{test_suite_main.yml => test_suite.yml} | 12 +-- .github/workflows/test_suite_release.yml | 24 ------ .gitignore | 3 - animate/adapt.py | 8 +- animate/cxx/metric2d.cxx | 2 +- animate/cxx/metric3d.cxx | 2 +- animate/cython/numbering.pyx | 59 ------------- animate/interpolation.py | 2 +- animate/math.py | 2 +- animate/metric.py | 20 ++--- animate/quality.py | 2 +- animate/recovery.py | 2 +- pyproject.toml | 34 ++------ setup.py | 83 ------------------- test/test_interpolation.py | 2 +- test/test_numbering.py | 47 ----------- test/test_quality.py | 2 +- test/test_recovery.py | 2 +- test/test_setup.py | 2 +- 19 files changed, 35 insertions(+), 275 deletions(-) rename .github/workflows/{test_suite_main.yml => test_suite.yml} (65%) delete mode 100644 .github/workflows/test_suite_release.yml delete mode 100644 animate/cython/numbering.pyx delete mode 100644 setup.py delete mode 100644 test/test_numbering.py diff --git a/.github/workflows/test_suite_main.yml b/.github/workflows/test_suite.yml similarity index 65% rename from .github/workflows/test_suite_main.yml rename to .github/workflows/test_suite.yml index 532108b4..d4f4ef5a 100644 --- a/.github/workflows/test_suite_main.yml +++ b/.github/workflows/test_suite.yml @@ -1,23 +1,23 @@ -name: 'Run Animate test suite (main)' +name: 'Run Animate Test Suite' on: - # Run test suite whenever the main branch is updated + # Run test suite whenever main is updated push: - branches: [main] + branches: + - main paths: - '**.py' - '**.cxx' - - '.github/workflows/test_suite_main.yml' + - '.github/workflows/test_suite.yml' - '.gitmodules' - 'pyproject.toml' # Run test suite whenever commits are pushed to an open PR pull_request: - branches: [main] paths: - '**.py' - '**.cxx' - - '.github/workflows/test_suite_main.yml' + - '.github/workflows/test_suite.yml' - '.gitmodules' - 'pyproject.toml' diff --git a/.github/workflows/test_suite_release.yml b/.github/workflows/test_suite_release.yml deleted file mode 100644 index 2b7e436e..00000000 --- a/.github/workflows/test_suite_release.yml +++ /dev/null @@ -1,24 +0,0 @@ -name: 'Run Animate test suite (release)' - -on: - # Run test suite whenever the release branch is updated - push: - branches: [release] - paths: - - '.github/workflows/test_suite_release.yml' - - # Run test suite whenever commits that change this workflow are pushed to an open PR - pull_request: - branches: [release] - paths: - - '.github/workflows/test_suite_release.yml' - - # Run test suite at 3AM on the first of the month - schedule: - - cron: '0 3 1 * *' - -jobs: - test_suite: - uses: mesh-adaptation/docs/.github/workflows/reusable_test_suite.yml@main - with: - docker-image: 'firedrake-parmmg:release' diff --git a/.gitignore b/.gitignore index 3e77c7bd..66b6cac1 100644 --- a/.gitignore +++ b/.gitignore @@ -10,12 +10,9 @@ __pycache__/ .coverage # File extensions -*.c -*.h5 *.jpg *.jpeg *.pdf *.png *.pvd -*.so *.vtu diff --git a/animate/adapt.py b/animate/adapt.py index 601164a9..7bd25c82 100644 --- a/animate/adapt.py +++ b/animate/adapt.py @@ -1,17 +1,17 @@ import abc import os -from functools import cached_property from shutil import rmtree import firedrake.checkpointing as fchk import firedrake.functionspace as ffs import firedrake.mesh as fmesh +import firedrake.utils as futils from firedrake import COMM_SELF, COMM_WORLD +from firedrake.cython.dmcommon import to_petsc_local_numbering from firedrake.petsc import PETSc from firedrake.projection import Projector from .checkpointing import get_checkpoint_dir, load_checkpoint, save_checkpoint -from .cython.numbering import to_petsc_local_numbering from .metric import RiemannianMetric __all__ = ["MetricBasedAdaptor", "adapt"] @@ -87,7 +87,7 @@ def __init__(self, mesh, metric, name=None, comm=None): self.metric = metric self.projectors = [] - @cached_property + @futils.cached_property @PETSc.Log.EventDecorator() def adapted_mesh(self): """ @@ -176,7 +176,7 @@ def adapt(mesh, *metrics, name=None, serialise=None, remove_checkpoints=True): """ nprocs = COMM_WORLD.size - dim = mesh.topological_dimension() + dim = mesh.topological_dimension if serialise is None: serialise = nprocs > 1 and dim != 3 elif not serialise and dim != 3: diff --git a/animate/cxx/metric2d.cxx b/animate/cxx/metric2d.cxx index 3e51c0aa..86d1b3e1 100644 --- a/animate/cxx/metric2d.cxx +++ b/animate/cxx/metric2d.cxx @@ -146,7 +146,7 @@ void intersect(double M_[4], const double * A_, const double * B_) { // Solve eigenvalue problem for triple product of inverse square root metric and the second metric SelfAdjointEigenSolver> eigensolver2(Sqi.transpose() * B * Sqi); Q = eigensolver2.eigenvectors(); - D = eigensolver2.eigenvalues().array().cwiseMax(1).matrix().asDiagonal(); + D = eigensolver2.eigenvalues().array().max(1).matrix().asDiagonal(); // Compute metric intersection M = Sq.transpose() * Q * D * Q.transpose() * Sq; diff --git a/animate/cxx/metric3d.cxx b/animate/cxx/metric3d.cxx index 8bc60eb2..78e0cc71 100644 --- a/animate/cxx/metric3d.cxx +++ b/animate/cxx/metric3d.cxx @@ -187,7 +187,7 @@ void intersect(double M_[9], const double * A_, const double * B_) { // Solve eigenvalue problem for triple product of inverse square root metric and the second metric SelfAdjointEigenSolver> eigensolver2(Sqi.transpose() * B * Sqi); Q = eigensolver2.eigenvectors(); - D = eigensolver2.eigenvalues().array().cwiseMax(1).matrix().asDiagonal(); + D = eigensolver2.eigenvalues().array().max(1).matrix().asDiagonal(); // Compute metric intersection M = Sq.transpose() * Q * D * Q.transpose() * Sq; diff --git a/animate/cython/numbering.pyx b/animate/cython/numbering.pyx deleted file mode 100644 index 3abaee2b..00000000 --- a/animate/cython/numbering.pyx +++ /dev/null @@ -1,59 +0,0 @@ -""" -Cython module for handling numberings of PETSc Vecs corresponding to Firedrake Functions. -""" -import numpy as np -from firedrake.petsc import PETSc -from firedrake.utils import IntType, ScalarType - -cimport numpy as np -cimport petsc4py.PETSc as PETSc - -from petsc4py.PETSc cimport CHKERR - - -cdef extern from "petsc.h": - ctypedef long PetscInt - ctypedef enum PetscErrorCode: - PETSC_SUCCESS - PETSC_ERR_LIB - - -cdef extern from "petscis.h" nogil: - PetscErrorCode PetscSectionGetDof(PETSc.PetscSection,PetscInt,PetscInt*) - PetscErrorCode PetscSectionGetOffset(PETSc.PetscSection,PetscInt,PetscInt*) - - -def to_petsc_local_numbering(PETSc.Vec vec, V): - """ - Reorder a PETSc Vec corresponding to a Firedrake Function w.r.t. - the PETSc natural numbering. - - :arg vec: the PETSc Vec to reorder; must be a global vector - :arg V: the FunctionSpace of the Function which the Vec comes from - :ret out: a copy of the Vec, ordered with the PETSc natural numbering - """ - cdef int dim, idx, start, end, p, d, k - cdef PetscInt dof, off - cdef PETSc.Vec out - cdef PETSc.Section section - cdef np.ndarray varray, oarray - - section = V.dm.getGlobalSection() - out = vec.duplicate() - varray = vec.array_r - oarray = out.array - dim = V.value_size - idx = 0 - start, end = vec.getOwnershipRange() - for p in range(*section.getChart()): - CHKERR(PetscSectionGetDof(section.sec, p, &dof)) - if dof > 0: - CHKERR(PetscSectionGetOffset(section.sec, p, &off)) - assert off >= 0 - off *= dim - for d in range(dof): - for k in range(dim): - oarray[idx] = varray[off + dim * d + k - start] - idx += 1 - assert idx == (end - start) - return out diff --git a/animate/interpolation.py b/animate/interpolation.py index 75423ef7..a7359390 100644 --- a/animate/interpolation.py +++ b/animate/interpolation.py @@ -324,7 +324,7 @@ def clement_interpolant(source, target_space=None, boundary=False): if rank not in (0, 1, 2): raise ValueError(f"Rank-{rank + 1} tensors are not supported.") mesh = Vs.mesh() - dim = mesh.topological_dimension() + dim = mesh.topological_dimension # Process target space Vt = target_space diff --git a/animate/math.py b/animate/math.py index 403723a1..4d492b15 100644 --- a/animate/math.py +++ b/animate/math.py @@ -73,7 +73,7 @@ def construct_basis(vector, normalise=True): if not isinstance(vector, ufl.core.expr.Expr): raise TypeError(f"Expected UFL Expr, not '{type(vector)}'.") as_vector = ufl.as_vector - dim = ufl.domain.extract_unique_domain(vector).topological_dimension() + dim = ufl.domain.extract_unique_domain(vector).topological_dimension if dim not in (2, 3): raise ValueError(f"Dimension {dim} not supported.") diff --git a/animate/metric.py b/animate/metric.py index e30d08db..2119fc6e 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -89,7 +89,7 @@ def __init__(self, function_space, *args, **kwargs): # Check that we have an appropriate tensor P1 function fs = self.function_space() mesh = fs.mesh() - tdim = mesh.topological_dimension() + tdim = mesh.topological_dimension if tdim not in (2, 3): raise ValueError(f"Riemannian metric should be 2D or 3D, not {tdim}D.") if isinstance(fs.dof_count, Iterable): @@ -303,7 +303,7 @@ def _set_plex_coordinates(self): consistent. """ entity_dofs = np.zeros(self._tdim + 1, dtype=np.int32) - entity_dofs[0] = self._mesh.geometric_dimension() + entity_dofs[0] = self._mesh.geometric_dimension coord_section = self._mesh.create_section(entity_dofs)[0] # NOTE: section doesn't have any fields, but PETSc assumes it to have one coord_dm = self._plex.getCoordinateDM() @@ -518,7 +518,7 @@ def interp(f): if not np.isclose(_a_max, 1.0) and _a_max < 1.0: raise ValueError(f"Encountered a_max value smaller than unity: {_a_max}.") - dim = mesh.topological_dimension() + dim = mesh.topological_dimension boundary_tag = self._variable_parameters.get("dm_plex_metric_boundary_tag") if boundary_tag is None: node_set = self.function_space().node_set @@ -734,7 +734,7 @@ def compute_eigendecomposition(self, reorder=False): mesh = V_ten.mesh() fe = (V_ten.ufl_element().family(), V_ten.ufl_element().degree()) V_vec = firedrake.VectorFunctionSpace(mesh, *fe) - dim = mesh.topological_dimension() + dim = mesh.topological_dimension evectors, evalues = firedrake.Function(V_ten), firedrake.Function(V_vec) if reorder: name = "get_reordered_eigendecomposition" @@ -788,7 +788,7 @@ def assemble_eigendecomposition(self, evectors, evalues): "Mismatching finite element space degrees:" f" {fe_ten.degree()} vs. {fe_vec.degree()}." ) - dim = V_ten.mesh().topological_dimension() + dim = V_ten.mesh().topological_dimension op2.par_loop( get_metric_kernel("set_eigendecomposition", dim), V_ten.node_set, @@ -848,7 +848,7 @@ def density_and_quotients(self, reorder=False): fs_ten = self.function_space() mesh = fs_ten.mesh() fe = (fs_ten.ufl_element().family(), fs_ten.ufl_element().degree()) - dim = mesh.topological_dimension() + dim = mesh.topological_dimension evectors, evalues = self.compute_eigendecomposition(reorder=reorder) # Extract density and quotients @@ -886,7 +886,7 @@ def compute_isotropic_metric( mesh = ufl.domain.extract_unique_domain(error_indicator) if mesh != self.function_space().mesh(): raise ValueError("Cannot use an error indicator from a different mesh.") - dim = mesh.topological_dimension() + dim = mesh.topological_dimension # Interpolate P0 indicators into P1 space if interpolant == "Clement": @@ -980,7 +980,7 @@ def compute_anisotropic_dwr_metric( mesh = ufl.domain.extract_unique_domain(error_indicator) if mesh != self.function_space().mesh(): raise ValueError("Cannot use an error indicator from a different mesh.") - dim = mesh.topological_dimension() + dim = mesh.topological_dimension if convergence_rate < 1.0: raise ValueError( f"Convergence rate must be at least one, not {convergence_rate}." @@ -1133,7 +1133,7 @@ def determine_metric_complexity(H_interior, H_boundary, target, p, **kwargs): :returns: unique solution of algebraic problem :rtype: :class:`float` """ - d = H_interior.function_space().mesh().topological_dimension() + d = H_interior.function_space().mesh().topological_dimension if d not in (2, 3): raise ValueError(f"Spatial dimension {d} not supported.") if np.isinf(p): @@ -1183,7 +1183,7 @@ def intersect_on_boundary(*metrics, boundary_tag="on_boundary"): n = len(metrics) assert n > 0, "Nothing to combine" fs = metrics[0].function_space() - dim = fs.mesh().topological_dimension() + dim = fs.mesh().topological_dimension if dim not in (2, 3): raise ValueError( f"Spatial dimension {dim} not supported. Must be either 2 or 3." diff --git a/animate/quality.py b/animate/quality.py index c05402d1..d3e060f4 100644 --- a/animate/quality.py +++ b/animate/quality.py @@ -60,7 +60,7 @@ def __init__(self, mesh, metric=None, python=False): self.mesh = mesh self.metric = metric self.python = python - self.dim = mesh.topological_dimension() + self.dim = mesh.topological_dimension self.coords = mesh.coordinates self.P0 = firedrake.FunctionSpace(mesh, "DG", 0) src_dir = os.path.join(os.path.dirname(__file__), "cxx") diff --git a/animate/recovery.py b/animate/recovery.py index 5b1e4577..e5ca1b5b 100644 --- a/animate/recovery.py +++ b/animate/recovery.py @@ -132,7 +132,7 @@ def recover_boundary_hessian(f, method="Clement", target_space=None, **kwargs): """ mesh = ufl.domain.extract_unique_domain(f) - d = mesh.topological_dimension() + d = mesh.topological_dimension assert d in (2, 3) # Apply Gram-Schmidt to get tangent vectors diff --git a/pyproject.toml b/pyproject.toml index 8f3e0d11..581e4e99 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,18 +1,9 @@ [build-system] -build-backend = "setuptools.build_meta" -requires = [ - "Cython", - "firedrake", - "numpy", - "petsc4py", - "petsctools", - "setuptools>=77.0.3", - "wheel", -] +requires = ["setuptools>=77.0.3"] [project] name = "animate" -version = "0.4" +version = "0.3" authors = [ {name = "Joseph G. Wallwork", email = "joe.wallwork@outlook.com"}, {name = "Davor Dundovic"}, @@ -52,25 +43,11 @@ Homepage = "https://mesh-adaptation.github.io" Documentation = "https://mesh-adaptation.github.io/animate/index.html" Repository = "https://github.com/mesh-adaptation/animate" -# Note: -# - Extension modules are generated dynamically at build time using the legacy -# helper (`setup.py`) which performs PETSc detection and populates Extension -# objects based on the local PETSc installation (via petsctools/petsc4py). -# - We intentionally do not declare ext-modules statically here so the build -# helper can compute include/library dirs and other flags based on an -# existing PETSc build. -# -# The build backend requirements above ensure that, during isolated builds, -# the build environment will have Cython, numpy, petsctools, petsc4py, etc. -# available before the build helper runs Cython to generate .c files. [tool.setuptools] packages = ["animate", "adapt_common"] [tool.setuptools.package-data] -animate = [ - "cxx/*.cxx", - "cython/*.pyx", -] +animate = ["cxx/*.cxx"] [tool.setuptools.package-dir] adapt_common = "adapt_common/adapt_common" @@ -86,7 +63,6 @@ select = [ "F", # Pyflakes "I", # isort ] - [tool.ruff.lint.per-file-ignores] "demos/*" = [ "E402", # module level import not at top of file @@ -96,6 +72,6 @@ select = [ [tool.pytest.ini_options] filterwarnings = [ - "ignore:`np.bool8` is a deprecated alias for `np.bool_`*:DeprecationWarning", - "ignore:unable to find git revision*:UserWarning", + "ignore:`np.bool8` is a deprecated alias for `np.bool_`*:DeprecationWarning", + "ignore:unable to find git revision*:UserWarning", ] diff --git a/setup.py b/setup.py deleted file mode 100644 index bfd6d59f..00000000 --- a/setup.py +++ /dev/null @@ -1,83 +0,0 @@ -"""Setup script for building Cython extensions for the `animate` package. - -The script performs PETSc detection and populates Extension objects based on the local -PETSc installation (via petsctools/petsc4py). -""" - -import os -from dataclasses import dataclass, field - -import numpy as np -import petsc4py -import petsctools -from Cython.Build import cythonize -from setuptools import Extension, find_packages, setup - - -@dataclass -class ExternalDependency: - """Courtesy of Firedrake. See - https://github.com/firedrakeproject/firedrake/blob/main/setup.py - """ - - include_dirs: list[str] = field(default_factory=list, init=True) - extra_compile_args: list[str] = field(default_factory=list, init=True) - libraries: list[str] = field(default_factory=list, init=True) - library_dirs: list[str] = field(default_factory=list, init=True) - extra_link_args: list[str] = field(default_factory=list, init=True) - runtime_library_dirs: list[str] = field(default_factory=list, init=True) - - def __add__(self, other): - combined = {} - for f in self.__dataclass_fields__.keys(): - combined[f] = getattr(self, f) + getattr(other, f) - return self.__class__(**combined) - - def keys(self): - return self.__dataclass_fields__.keys() - - def __getitem__(self, key): - try: - return getattr(self, key) - except AttributeError as attr_err: - raise KeyError(f"Key {key} not present") from attr_err - - -def extensions(): - """Returns a list of Cython extensions to be compiled.""" - - # Define external dependencies - mpi_ = ExternalDependency( - extra_compile_args=petsctools.get_petscvariables()["MPICC_SHOW"].split()[1:], - ) - numpy_ = ExternalDependency(include_dirs=[np.get_include()]) - petsc_dir = petsctools.get_petsc_dir() - petsc_dirs = [petsc_dir, os.path.join(petsc_dir, petsctools.get_petsc_arch())] - petsc_includes = [petsc4py.get_include()] + [ - os.path.join(d, "include") for d in petsc_dirs - ] - petsc_ = ExternalDependency( - libraries=["petsc"], - include_dirs=petsc_includes, - library_dirs=[os.path.join(petsc_dirs[-1], "lib")], - runtime_library_dirs=[os.path.join(petsc_dirs[-1], "lib")], - ) - - # Define Cython extensions - cython_list = [ - Extension( - name="animate.cython.numbering", - language="c", - sources=[os.path.join("animate", "cython", "numbering.pyx")], - **(mpi_ + petsc_ + numpy_), - ) - ] - return cythonize(cython_list) - - -if __name__ == "__main__": - # Run the setup function to build the Cython extensions - setup( - packages=find_packages(), - ext_modules=extensions(), - ) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 67ad6d62..a02166d9 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -66,7 +66,7 @@ def get_space(self, rank, family, degree): elif rank == 1: return VectorFunctionSpace(self.mesh, family, degree) else: - shape = tuple(rank * [self.mesh.topological_dimension()]) + shape = tuple(rank * [self.mesh.topological_dimension]) return TensorFunctionSpace(self.mesh, family, degree, shape=shape) def analytic(self, rank): diff --git a/test/test_numbering.py b/test/test_numbering.py deleted file mode 100644 index 9d859c30..00000000 --- a/test/test_numbering.py +++ /dev/null @@ -1,47 +0,0 @@ -"""Module containing unit tests for the cython.numbering module.""" - -import firedrake as fd -import numpy as np -import pytest -from animate.cython.numbering import to_petsc_local_numbering - - -@pytest.fixture -def mesh(): - """Fixture to create a uniform unit square mesh.""" - return fd.UnitSquareMesh(4, 4) - - -@pytest.fixture(params=[0, 1, 2]) -def rank(request): - """Fixture specifying function space rank.""" - return request.param - - -@pytest.fixture -def function_space(mesh, rank): - """Fixture to create a function space on a given mesh with a given rank.""" - try: - return { - 0: fd.FunctionSpace, - 1: fd.VectorFunctionSpace, - 2: fd.TensorFunctionSpace, - }[rank](mesh, "CG", 1) - except KeyError as key_err: - raise ValueError(f"Rank {rank} not considered.") from key_err - - -def test_is_permutation(function_space): - """Test that to_petsc_local_numbering is indeed a permutation.""" - f = fd.Function(function_space) - - # Fill the vector with arbitrary values - with f.dat.vec as vec: - vec.array[:] = np.arange(vec.size) - - # Verify that reordering according to the PETSc numbering is a permutation - with f.dat.vec_ro as vec: - reordered_vec = to_petsc_local_numbering(vec, function_space) - assert reordered_vec is not None - assert reordered_vec.size == vec.size - assert sorted(reordered_vec.array) == sorted(vec.array) diff --git a/test/test_quality.py b/test/test_quality.py index 703b94b3..133163c9 100644 --- a/test/test_quality.py +++ b/test/test_quality.py @@ -37,7 +37,7 @@ class TestQuality(unittest.TestCase): """ def quality(self, name, mesh, **kwargs): - dim = mesh.topological_dimension() + dim = mesh.topological_dimension if name == "metric": P1_ten = TensorFunctionSpace(mesh, "CG", 1) M = Function(P1_ten).interpolate(ufl.Identity(dim)) diff --git a/test/test_recovery.py b/test/test_recovery.py index 39cd1a66..1e7ebe4e 100644 --- a/test/test_recovery.py +++ b/test/test_recovery.py @@ -120,7 +120,7 @@ def metric(mesh): @staticmethod def relative_error(approx, ignore_boundary=False): mesh = approx.function_space().mesh() - dim = mesh.topological_dimension() + dim = mesh.topological_dimension P1_ten = TensorFunctionSpace(mesh, "CG", 1) identity = Function(P1_ten).interpolate(ufl.Identity(dim)) diff --git a/test/test_setup.py b/test/test_setup.py index 29364fea..42589aca 100644 --- a/test/test_setup.py +++ b/test/test_setup.py @@ -49,7 +49,7 @@ def uniform_metric(mesh, a=100.0, metric_parameters=None): else: function_space = mesh mesh = function_space.mesh() - dim = mesh.topological_dimension() + dim = mesh.topological_dimension metric = RiemannianMetric(function_space) metric.interpolate(a * ufl.Identity(dim)) metric.set_parameters(metric_parameters)