DirichletBC control

CITATION.cff

@@ -20,5 +20,5 @@ keywords:
   - Adjoint
   - Optimization
 license: MIT
-version: v0.1.0
-date-released: '2023-04-04'
+version: v0.3.0
+date-released: '2026-08-06'

README.md

@@ -52,6 +52,15 @@ cd dolfinx-adjoint
 python3 -m pip install -e ".[all]"
 ```
 
+The `docs` extra builds the demos and API docs, but not [Moola](https://github.com/funsim/moola)
+(used by two of the three demos): PyPI rejects packages that declare a direct git-URL
+dependency, so it can't be listed in `pyproject.toml`. Install it separately first if you
+want to build the docs or run those demos locally:
+
+```bash
+python3 -m pip install git+https://github.com/funsim/moola.git
+```
+
 
 ## Quick Start
 

demos/demo_nonmatching_grids.py

@@ -6,7 +6,7 @@
 #
 # Subject to the state equation on a fine grid (\Omega_s) and a control on a coarse grid (\Omega_c):
 #
-# $$ - \Delta y = I_{\Omega_c \to \Omega_s}(u + u^3) \quad \text{in } \Omega_s $$
+# $$ - \Delta y = I_{\Omega_c \to \Omega_s}(u + 0.1 u^3) \quad \text{in } \Omega_s $$
 
 from mpi4py import MPI
 

pyproject.toml

@@ -10,7 +10,7 @@ license = "MIT"
 license-files = ["LICENSE"]
 readme = "README.md"
 dependencies = [
-      "fenics-dolfinx>=0.10.0",
+      "fenics-dolfinx>=0.11.0",
       "pyadjoint-ad>=2025.10.0",
       "typing_extensions; python_version < '3.11'",
       "packaging>=24.2",
@@ -25,7 +25,10 @@ dev = ["pdbpp", "ipython", "mypy", "ruff"]
 docs = [
       "jupyter-book<2.0",
       "jupytext",
-      # "moola@git+https://github.com/funsim/moola.git",
+      # moola (needed by two of the three demos) is deliberately not listed here: PyPI
+      # #     mesh = form.ufl_domain()
+        #     space = c1._ad_function_space(mesh)rejects packages that declare a direct git-URL dependency. Install it manually
+      # first (see README.md's Development Install section, and .github/workflows/build_docs.yml).
       "pandas",
       "pyvista[all]>0.45",
       "networkx",
@@ -102,3 +105,8 @@ current_version = "0.3.0"
 filename = "pyproject.toml"
 search = 'version = "{current_version}"'
 replace = 'version = "{new_version}"'
+
+[[tool.bumpversion.files]]
+filename = "CITATION.cff"
+search = 'version: v{current_version}'
+replace = 'version: v{new_version}'

src/dolfinx_adjoint/__init__.py

@@ -1,4 +1,4 @@
-"""Top-level package for dxa."""
+"""Top-level package for dolfinx_adjoint."""
 
 from importlib.metadata import metadata
 

src/dolfinx_adjoint/assembly.py

@@ -21,8 +21,8 @@ def assemble_scalar(form: ufl.Form, **kwargs):
             Includes ``"ad_block_tag"`` to tag the block in the adjoint tape,
             ``"annotate"`` to control whether the assembly is annotated in the adjoint tape,
             ``"jit_options"`` for JIT compilation options,
-            and ``"form_compiler_options"`` for form compiler options and ``"entity_map"`` for assembling with Arguments
-            and coefficients form meshes that has some relation.
+            ``"form_compiler_options"`` for form compiler options, and ``"entity_maps"`` for
+            assembling with Arguments and coefficients form meshes that has some relation.
     """
     ad_block_tag = kwargs.pop("ad_block_tag", None)
 
@@ -56,7 +56,7 @@ def assemble_scalar(form: ufl.Form, **kwargs):
 def error_norm(
     u_ex: ufl.core.expr.Expr,
     u: ufl.core.expr.Expr,
-    norm_type=typing.Literal["L2", "H1"],
+    norm_type: typing.Literal["L2", "H1"] = "L2",
     jit_options: dict | None = None,
     form_compiler_options: dict | None = None,
     entity_map: dict[dolfinx.mesh.Mesh, npt.NDArray[numpy.int32]] | None = None,

src/dolfinx_adjoint/blocks/__init__.py

@@ -1,7 +1,17 @@
 from .assembly import AssembleBlock
+from .dirichletbc import DirichletBCBlock
 from .function_assigner import FunctionAssignBlock
+from .interpolation import ExprInterpolationBlock, InterpolationBlock
+from .nonmatching_interpolation import NonmatchingInterpolationBlock
+from .solvers import LinearProblemBlock, NonlinearProblemBlock
 
 __all__ = [
     "AssembleBlock",
+    "DirichletBCBlock",
+    "ExprInterpolationBlock",
     "FunctionAssignBlock",
+    "InterpolationBlock",
+    "LinearProblemBlock",
+    "NonlinearProblemBlock",
+    "NonmatchingInterpolationBlock",
 ]

src/dolfinx_adjoint/blocks/_vector.py

@@ -36,6 +36,7 @@ def _vector(
         map: Index map the describes the size and distribution of the
             vector.
         bs: Block size.
+        function_space: The function space the vector is associated with.
         dtype: The scalar type.
 
     Returns:
@@ -75,6 +76,8 @@ def _create_vector(L: dolfinx.fem.Form, space: dolfinx.fem.FunctionSpace) -> _Sp
 
     Args:
         L: A linear form.
+        space: The function space ``L``'s (single) argument lives on -- must match
+            ``L.function_spaces[0]``.
 
     Returns:
         A vector that the form can be assembled into.

src/dolfinx_adjoint/blocks/assembly.py

@@ -13,19 +13,22 @@
 def assemble_compiled_form(
     form: dolfinx.fem.Form, tensor: typing.Union[dolfinx.la.Vector, _SpecialVector | float] | None = None
 ) -> typing.Union[dolfinx.la.Vector, _SpecialVector, float]:
-    """Assemble a compiled form and optionally apply Dirichlet boundary condition.
+    """Assemble a compiled form into ``tensor`` (or return a new scalar).
 
     Args:
         form: Compiled form to assemble.
-        tensor: Optional vector to which the assembled form will be added.
+        tensor: For a rank-1 form, the vector to accumulate the assembled contribution
+            into, while it is unused for a rank-0 form.
     Returns:
-        tensor: The assembled vector, which is either the input tensor or a new vector
-            created from the form's function space(s).
+        For a rank-1 form, ``tensor`` itself (mutated in place). For a rank-0 form, the
+        assembled scalar as a Python ``float``.
     Raises:
         NotImplementedError: If the form's rank is not 0 or 1.
     """
 
     if form.rank == 1:
+        if tensor is None:
+            raise ValueError("tensor must be provided for rank-1 forms.")
         assert isinstance(tensor, dolfinx.la.Vector)
         dolfinx.fem.assemble._assemble_vector_array(tensor.array, form)
         tensor.scatter_reverse(dolfinx.la.InsertMode.add)
@@ -34,9 +37,9 @@ def assemble_compiled_form(
         local_val = dolfinx.fem.assemble_scalar(form)
         comm = form.mesh.comm
         tensor = comm.allreduce(local_val, op=MPI.SUM)
-
     else:
         raise NotImplementedError("Only 1-form assembly is currently supported.")
+    assert tensor is not None
     return tensor
 
 
@@ -271,17 +274,19 @@ def evaluate_hessian_component(
         c1_rep = block_variable.saved_output
 
         if isinstance(c1, dolfinx.fem.Constant):
-            mesh = form.ufl_domain()
-            space = c1._ad_function_space(mesh)
-        elif isinstance(c1, dolfinx.fem.Function):
+            raise RuntimeError(
+                "All constants should have been replaced with real space coefficients before this point."
+            )
+        if isinstance(c1, dolfinx.fem.Function):
             space = c1.function_space
-        elif isinstance(c1, dolfinx.mesh.Mesh):
-            c1_rep = ufl.SpatialCoordinate(c1)
-            space = c1._ad_function_space()
+        # TODO: Add support for shape optimization
+        # elif isinstance(c1, dolfinx.mesh.Mesh):
+        #     c1_rep = ufl.SpatialCoordinate(c1)
+        #     space = c1._ad_function_space()
         else:
             return None
         hessian_outputs, dform = self.compute_action_adjoint(hessian_input, arity_form, form, c1_rep, space)
-        ddform = 0
+        ddform = 0.0
         for other_idx, bv in relevant_dependencies:
             c2_rep = bv.saved_output
             tlm_input = bv.tlm_value
@@ -298,7 +303,6 @@ def evaluate_hessian_component(
             ddform = ufl.algorithms.expand_derivatives(ddform)
 
         if not ddform.empty():
-            # FIXME: COmpare ddform with legacy dolfin_adjoitn here, as this is DG-0, while hessian is in DG-0
             adj_action = self.compute_action_adjoint(adj_input, arity_form, dform=ddform)[0]
             try:
                 hessian_outputs += adj_action

src/dolfinx_adjoint/blocks/dirichletbc.py

@@ -1,9 +1,75 @@
 import dolfinx
 import numpy as np
 import numpy.typing as npt
+from packaging.version import Version
 from pyadjoint.block import Block
 
 
+def build_cpp_bc_and_kwargs(
+    g: dolfinx.fem.Function, dofs: npt.NDArray[np.int32], V: dolfinx.fem.FunctionSpace | None
+) -> tuple[
+    dolfinx.cpp.fem.DirichletBC_float32
+    | dolfinx.cpp.fem.DirichletBC_float64
+    | dolfinx.cpp.fem.DirichletBC_complex64
+    | dolfinx.cpp.fem.DirichletBC_complex128,
+    dict,
+]:
+    """Build the cpp Dirichlet bc object and the kwargs a Python
+    :py:class:`dolfinx.fem.DirichletBC` wrapper needs, decoupling the *Python-level*
+    constrained space ``V`` from the cpp-level construction.
+
+    Shared by :py:class:`~dolfinx_adjoint.types.dirichletbc.DirichletBC` (the tape-tracked
+    constructor) and :py:meth:`DirichletBCBlock.evaluate_tlm_component` (a plain,
+    untracked bc built fresh each TLM call).
+
+    The cpp constructor is always called the V-less way regardless of whether the caller
+    passed ``V``: the 3-arg cpp overload for a `Function`-valued ``g`` requires ``dofs`` to
+    be a *paired* ``(dofs_in_V, dofs_in_g_space)`` sequence -- the mechanism for
+    constraining a sub-space with a value on its collapsed counterpart, not what a flat
+    ``dofs`` array plus a broadcast-style value (e.g. a
+    :py:class:`~dolfinx_adjoint.Constant`, itself a `Function` on a single-dof real space)
+    needs. The V-less overload already broadcasts a real-space `Function`'s single value
+    across an arbitrary flat ``dofs`` array correctly. ``V`` is instead threaded through
+    purely at the *Python* level, which :py:class:`dolfinx.fem.DirichletBC` keeps entirely
+    independent of the cpp object: ``bc.function_space`` is whatever ``V`` is passed to the
+    Python wrapper, never introspected from the cpp bc. See dolfinx-adjoint-knowledge's
+    scratch/boundary-control/issues/01 for the bug this avoids.
+    """
+    dtype = g.dtype
+    V_used = V if V is not None else g.function_space
+    cpp_bc: (
+        dolfinx.cpp.fem.DirichletBC_float32
+        | dolfinx.cpp.fem.DirichletBC_float64
+        | dolfinx.cpp.fem.DirichletBC_complex64
+        | dolfinx.cpp.fem.DirichletBC_complex128
+    )
+    # cpp_bc is constructed inside each branch, immediately after its own isinstance
+    # assert, rather than deferred to one call after the if/elif chain -- mypy can only
+    # narrow g._cpp_object's union type within the branch the assert itself is in.
+    if np.issubdtype(dtype, np.float32):
+        assert isinstance(g._cpp_object, dolfinx.cpp.fem.Function_float32)
+        cpp_bc = dolfinx.cpp.fem.DirichletBC_float32(g._cpp_object, dofs)
+    elif np.issubdtype(dtype, np.float64):
+        assert isinstance(g._cpp_object, dolfinx.cpp.fem.Function_float64)
+        cpp_bc = dolfinx.cpp.fem.DirichletBC_float64(g._cpp_object, dofs)
+    elif np.issubdtype(dtype, np.complex64):
+        assert isinstance(g._cpp_object, dolfinx.cpp.fem.Function_complex64)
+        cpp_bc = dolfinx.cpp.fem.DirichletBC_complex64(g._cpp_object, dofs)
+    elif np.issubdtype(dtype, np.complex128):
+        assert isinstance(g._cpp_object, dolfinx.cpp.fem.Function_complex128)
+        cpp_bc = dolfinx.cpp.fem.DirichletBC_complex128(g._cpp_object, dofs)
+    else:
+        raise NotImplementedError(f"Type {dtype} not supported.")
+
+    bc_kwargs: dict = {}
+    # If dolfinx-version is 0.12 we need to pass the following
+    # due to https://github.com/FEniCS/dolfinx/pull/4342/
+    if Version(dolfinx.__version__).minor > 11:
+        bc_kwargs["V"] = V_used
+        bc_kwargs["g"] = g
+    return cpp_bc, bc_kwargs
+
+
 class DirichletBCBlock(Block):
     """A block representing a DirichletBC in the adjoint framework.
 
@@ -39,4 +105,65 @@ def prepare_recompute_component(self, inputs, relevant_outputs):
         return inputs[0] if inputs else None
 
     def recompute_component(self, inputs, block_variable, idx, prepared):
-        return block_variable.saved_output
+        """Return the (aliased) bc, having first resynced its live ``g.x.array`` from the
+        value's own checkpoint at this tape position.
+
+        ``DirichletBC._ad_create_checkpoint``/``_ad_restore_at_checkpoint``
+        (types/dirichletbc.py) both ``return self`` -- the bc's own "checkpoint" aliases
+        the *live* bc object rather than snapshotting a value -- so nothing else writes a
+        replayed/perturbed value back into ``bc.g``'s array. ``prepared`` (this block's
+        single dependency's own, correctly-checkpointed value, from
+        ``prepare_recompute_component``) is exactly that value: writing it into ``bc.g``
+        here, at the position in the tape this block itself occupies (always *before* any
+        solve block that consumes ``bc``, since the bc must be constructed first), is what
+        makes a later solve block see the right value regardless of which tape position is
+        being replayed. See dolfinx-adjoint-knowledge's
+        scratch/boundary-control/issues/02 for what breaks without this.
+        """
+        bc = block_variable.saved_output
+        if prepared is not None:
+            bc.g.x.array[:] = prepared.x.array[:]
+            bc.g.x.scatter_forward()
+        return bc
+
+    def evaluate_tlm_component(self, inputs, tlm_inputs, block_variable, idx, prepared=None):
+        """Return this bc's own tangent-linear value: itself a (plain, untracked)
+        DirichletBC, with its value replaced by the perturbation direction.
+
+        A bc perturbation is *not* an ordinary right-hand-side contribution -- it enters
+        the tangent-linear solve as an inhomogeneous condition (`u_dot = g_dot` on this
+        bc's dofs, see ``HomogeneousBCLinearProblem.tlm_bcs``/``solve()``), consumed by
+        ``_ProblemBlockBase.prepare_evaluate_tlm``. Built plain (``dolfinx.fem.dirichletbc``,
+        not the overloaded ``dolfinx_adjoint`` one) so it is never itself tape-recorded --
+        it exists only for this one TLM evaluation, not as a new control.
+        """
+        tlm_input = tlm_inputs[0]
+        if tlm_input is None:
+            return None
+        cpp_bc, bc_kwargs = build_cpp_bc_and_kwargs(tlm_input, self._dofs, self._V)
+        return dolfinx.fem.DirichletBC(cpp_bc, **bc_kwargs)
+
+    def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
+        """Return this bc's contribution to the adjoint action.
+
+        ``adj_inputs[0]`` is the boundary reaction the consuming solve block(s) already
+        computed (see ``_ProblemBlockBase._mask_reaction_to_bc``/``prepare_evaluate_adj``),
+        already living on this bc's own constrained space -- no reduction needed here:
+        ``types.dirichletbc._pack_bc_value`` always packs the bc's value into a Function on
+        exactly that space before this block is ever created, whatever the original value
+        was (a plain `Function`, a broadcasting `Constant`, or a general expression), so
+        this block's single dependency and the masked reaction always agree.
+        """
+        return adj_inputs[0]
+
+    def evaluate_hessian_component(
+        self, inputs, hessian_inputs, adj_inputs, block_variable, idx, relevant_dependencies, prepared=None
+    ):
+        """Return this bc's contribution to the Hessian action.
+
+        Same pass-through as ``evaluate_adj_component``, applied to the second-order
+        boundary reaction (``prepare_evaluate_hessian``'s ``_adj_sol2_bdy``) instead of the
+        first-order one -- the *entire* Hessian-action contribution for a Dirichlet bc
+        control (see dolfinx-adjoint-knowledge's scratch/boundary-control/spec.md for why).
+        """
+        return hessian_inputs[0]

src/dolfinx_adjoint/blocks/interpolation.py

@@ -436,6 +436,16 @@ def prepare_evaluate_hessian(self, inputs, hessian_inputs, adj_inputs, relevant_
                     V_in = target_dep_i.function_space
                     du = ufl.TrialFunction(V_in)
                     d2E = ufl.derivative(dE_total, target_dep_i, du)
+                    # ufl.derivative returns a lazy, unexpanded CoefficientDerivative node
+                    # that formally references `du` regardless of whether the expanded
+                    # expression actually depends on it (e.g. `dE_total` linear in
+                    # target_dep_i, as for a bare-coefficient expr -- its second
+                    # derivative is identically zero, but the *unexpanded* node still
+                    # reports one argument, previously causing a spurious H_op to be
+                    # compiled from a mesh-less zero expression). Expand derivatives
+                    # first so the argument count (and isinstance-zero check) reflect
+                    # the true, simplified expression.
+                    d2E = ufl.algorithms.apply_derivatives.apply_derivatives(d2E)
 
                     if not isinstance(d2E, (int, float)):
                         args = ufl.algorithms.extract_arguments(d2E)