Add P argument to solve function

.github/workflows/build.yml

@@ -23,9 +23,8 @@ env:
 jobs:
   build:
     name: "Build Firedrake"
-    # The type of runner that the job will run on
+    # Run on our self-hosted machines
     runs-on: self-hosted
-    # The docker container to use.
     container:
       image: firedrakeproject/firedrake-env:latest
     strategy:
@@ -49,7 +48,7 @@ jobs:
       COMPLEX: ${{ matrix.complex }}
       RDMAV_FORK_SAFE: 1
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Cleanup
         if: ${{ always() }}
         run: |

.github/workflows/docker_reuse.yml

@@ -29,7 +29,7 @@ jobs:
     runs-on: self-hosted
     steps:
       - name: Check out the repo
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
       - name: Log in to Docker Hub
         uses: docker/login-action@v2
         with:

.github/workflows/docs.yml

@@ -18,18 +18,19 @@ concurrency:
 jobs:
   build_docs:
     name: Run doc build
-    # The type of runner that the job will run on
+    # Run on the Github hosted runner
     runs-on: ubuntu-latest
-    # The docker container to use.
     container:
       image: firedrakeproject/firedrake-docdeps:latest
+      # Github hosted runners require running as root user:
+      # https://docs.github.com/en/actions/using-github-hosted-runners/about-github-hosted-runners/about-github-hosted-runners#docker-container-filesystem
       options: --user root
       volumes:
         - ${{ github.workspace }}:/home/firedrake/output
     # Steps represent a sequence of tasks that will be executed as
     # part of the jobs
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Install checkedout Firedrake
         run: |
           . /home/firedrake/firedrake/bin/activate
@@ -57,7 +58,7 @@ jobs:
           cd docs
           cp build/latex/Firedrake.pdf build/html/_static/manual.pdf
       - name: Upload artifact
-        uses: actions/upload-pages-artifact@v1
+        uses: actions/upload-pages-artifact@v3
         with:
           name: github-pages
           path: /__w/firedrake/firedrake/docs/build/html
@@ -76,4 +77,4 @@ jobs:
     steps:
       - name: Deploy to GitHub Pages
         id: deployment
-        uses: actions/deploy-pages@v2
+        uses: actions/deploy-pages@v4

.github/workflows/lint.yml

@@ -12,7 +12,7 @@ jobs:
     name: "Run linter"
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Setup python
         uses: actions/setup-python@v4
         with:
@@ -40,7 +40,7 @@ jobs:
     # #example-error-annotation-on-github-actions
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Check workflow files
         uses: docker://rhysd/actionlint:latest
         with:
@@ -53,7 +53,7 @@ jobs:
       env:
         HADOLINT_IGNORE: "DL3005,DL3007,DL3008,DL3015,DL3059"
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Lint dockerfiles inside hadolint container
         run: |
           for DOCKERFILE in docker/Dockerfile.*; \

.github/workflows/zenodo-canary.yml

@@ -13,7 +13,7 @@ jobs:
     name: "Run zenodo canary"
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
         with:
           fetch-depth: 1
       - name: Setup Python
@@ -27,7 +27,7 @@ jobs:
         run: |
           python scripts/firedrake-install --test-doi-resolution
       - name: Upload log
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v4
         if: failure()
         with:
           name: "zenodo-canary"

firedrake/adjoint_utils/function.py

@@ -275,8 +275,13 @@ def _define_riesz_map_form(self, riesz_representation, V):
 
     @no_annotations
     def _ad_convert_type(self, value, options=None):
-        # `_ad_convert_type` is not annoated unlike to `_ad_convert_riesz`
-        return self._ad_convert_riesz(value, options=options)
+        # `_ad_convert_type` is not annotated, unlike `_ad_convert_riesz`
+        options = {} if options is None else options
+        riesz_representation = options.get("riesz_representation", "L2")
+        if riesz_representation is None:
+            return value
+        else:
+            return self._ad_convert_riesz(value, options=options)
 
     def _ad_restore_at_checkpoint(self, checkpoint):
         if isinstance(checkpoint, CheckpointBase):

firedrake/functionspaceimpl.py

@@ -50,8 +50,8 @@ def check_element(element, top=True):
         If the element is illegal.
     """
     if element.cell.cellname() == "hexahedron" and \
-       element.family() not in ["Q", "DQ"]:
-        raise NotImplementedError("Currently can only use 'Q' and/or 'DQ' elements on hexahedral meshes, not", element.family())
+       element.family() not in ["Q", "DQ", "Real"]:
+        raise NotImplementedError("Currently can only use 'Q', 'DQ', and/or 'Real' elements on hexahedral meshes, not", element.family())
     if type(element) in (finat.ufl.BrokenElement, finat.ufl.RestrictedElement,
                          finat.ufl.HDivElement, finat.ufl.HCurlElement):
         inner = (element._element, )

firedrake/preconditioners/fdm.py

@@ -20,6 +20,7 @@
 from ufl.algorithms.expand_indices import expand_indices
 from tsfc.finatinterface import create_element
 from pyop2.compilation import load
+from pyop2.mpi import COMM_SELF
 from pyop2.sparsity import get_preallocation
 from pyop2.utils import get_petsc_dir, as_tuple
 from pyop2 import op2
@@ -226,9 +227,9 @@ def allocate_matrix(self, Amat, V, J, bcs, fcp, pmat_type, use_static_condensati
         self.embedding_element = ebig
 
         if Vbig.value_size == 1:
-            self.fises = PETSc.IS().createGeneral(fdofs, comm=PETSc.COMM_SELF)
+            self.fises = PETSc.IS().createGeneral(fdofs, comm=COMM_SELF)
         else:
-            self.fises = PETSc.IS().createBlock(Vbig.value_size, fdofs, comm=PETSc.COMM_SELF)
+            self.fises = PETSc.IS().createBlock(Vbig.value_size, fdofs, comm=COMM_SELF)
 
         # Create data structures needed for assembly
         self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V}
@@ -626,7 +627,7 @@ def assemble_reference_tensor(self, V, transpose=False, sort_interior=False):
                 q1 = sorted(get_base_elements(Q1), key=lambda e: e.formdegree)[-1]
 
             # Interpolate V * d(V) -> space(beta) * space(alpha)
-            comm = PETSc.COMM_SELF
+            comm = COMM_SELF
             zero = PETSc.Mat()
             A00 = petsc_sparse(evaluate_dual(e0, q0), comm=comm) if e0 and q0 else zero
             A11 = petsc_sparse(evaluate_dual(e1, q1), comm=comm) if e1 else zero
@@ -661,7 +662,7 @@ def _element_mass_matrix(self):
         if shape[2] > 1:
             ai *= shape[2]
             data = numpy.tile(numpy.eye(shape[2], dtype=data.dtype), shape[:1] + (1,)*(len(shape)-1))
-        return PETSc.Mat().createAIJ((nrows, nrows), csr=(ai, aj, data), comm=PETSc.COMM_SELF)
+        return PETSc.Mat().createAIJ((nrows, nrows), csr=(ai, aj, data), comm=COMM_SELF)
 
     @PETSc.Log.EventDecorator("FDMSetValues")
     def set_values(self, A, Vrow, Vcol, addv, mat_type="aij"):
@@ -1594,9 +1595,9 @@ def tabulate_exterior_derivative(Vc, Vf, cbcs=[], fbcs=[], comm=None):
 
     tdim = Vc.mesh().topological_dimension()
     zero = PETSc.Mat()
-    A00 = petsc_sparse(evaluate_dual(c0, f0), comm=PETSc.COMM_SELF) if f0 else zero
-    A11 = petsc_sparse(evaluate_dual(c1, f1), comm=PETSc.COMM_SELF) if c1 else zero
-    A10 = petsc_sparse(evaluate_dual(c0, f1, "grad"), comm=PETSc.COMM_SELF)
+    A00 = petsc_sparse(evaluate_dual(c0, f0), comm=COMM_SELF) if f0 else zero
+    A11 = petsc_sparse(evaluate_dual(c1, f1), comm=COMM_SELF) if c1 else zero
+    A10 = petsc_sparse(evaluate_dual(c0, f1, "grad"), comm=COMM_SELF)
     Dhat = block_mat(diff_blocks(tdim, ec.formdegree, A00, A11, A10), destroy_blocks=True)
     A00.destroy()
     A11.destroy()
@@ -1822,7 +1823,7 @@ def assemble_reference_tensor(self, V):
             try:
                 rtensor = cache[key]
             except KeyError:
-                rtensor = cache.setdefault(key, fdm_setup_ipdg(e, eta, comm=PETSc.COMM_SELF))
+                rtensor = cache.setdefault(key, fdm_setup_ipdg(e, eta, comm=COMM_SELF))
             Afdm[:0], Dfdm[:0], bdof[:0] = tuple(zip(rtensor))
             if not is_dg and e.degree() == degree:
                 # do not apply SIPG along continuous directions
@@ -1849,7 +1850,7 @@ def set_values(self, A, Vrow, Vcol, addv, mat_type="aij"):
             to determine the nonzeros on the upper triangual part of an ``'sbaij'`` matrix.
         """
         triu = A.getType() == "preallocator" and mat_type.endswith("sbaij")
-        set_submat = SparseAssembler.setSubMatCSR(PETSc.COMM_SELF, triu=triu)
+        set_submat = SparseAssembler.setSubMatCSR(COMM_SELF, triu=triu)
         update_A = lambda A, Ae, rindices: set_submat(A, Ae, rindices, rindices, addv)
         condense_element_mat = lambda x: x
 
@@ -1907,7 +1908,7 @@ def cell_to_global(lgmap, cell_to_local, cell_index, result=None):
             for e in range(nel):
                 # Ae = Be kron Bq[e]
                 adata = numpy.sum(Bq.dat.data_ro[index_coef(e)], axis=0)
-                Ae = PETSc.Mat().createAIJWithArrays(bshape, (aptr, aidx, adata), comm=PETSc.COMM_SELF)
+                Ae = PETSc.Mat().createAIJWithArrays(bshape, (aptr, aidx, adata), comm=COMM_SELF)
                 Ae = Be.kron(Ae)
                 rindices = get_rindices(e, result=rindices)
                 update_A(A, Ae, rindices)
@@ -2244,7 +2245,7 @@ def fdm_setup_ipdg(fdm_element, eta, comm=None):
 
     :arg fdm_element: a :class:`FIAT.FDMElement`
     :arg eta: penalty coefficient as a `float`
-    :arg comm: a :class:`PETSc.Comm`
+    :arg comm: an mpi4py communicator
 
     :returns: 3-tuple of:
         Afdm: a list of :class:`PETSc.Mats` with the sparse interval matrices

firedrake/solving.py

@@ -47,11 +47,19 @@ def solve(*args, **kwargs):
 
         solve(A, x, b, bcs=bcs, solver_parameters={...})
 
-    where `A` is a :class:`.Matrix` and `x` and `b` are :class:`.Function`\s.
-    If present, `bcs` should be a list of :class:`.DirichletBC`\s and
+    where ``A`` is a :class:`.Matrix` and ``x`` and ``b`` are :class:`.Function`\s.
+    If present, ``bcs`` should be a list of :class:`.DirichletBC`\s and
     :class:`.EquationBC`\s specifying, respectively, the strong boundary conditions
     to apply and PDEs to solve on the boundaries.
-    For the format of `solver_parameters` see below.
+    For the format of ``solver_parameters`` see below.
+
+    Optionally, an argument ``P`` of type :class:`~.MatrixBase` can be passed to
+    construct any preconditioner from; if none is supplied ``A`` is used to
+    construct the preconditioner.
+
+    .. code-block:: python3
+
+        solve(A, x, b, P=P, bcs=bcs, solver_parameters={...})
 
     *2. Solving linear variational problems*
 
@@ -192,6 +200,9 @@ def _la_solve(A, x, b, **kwargs):
     :arg A: the assembled bilinear form, a :class:`.Matrix`.
     :arg x: the :class:`.Function` to write the solution into.
     :arg b: the :class:`.Function` defining the right hand side values.
+    :kwarg P: an optional :class:`~.MatrixBase` to construct any
+         preconditioner from; if none is supplied ``A`` is
+         used to construct the preconditioner.
     :kwarg solver_parameters: optional solver parameters.
     :kwarg nullspace: an optional :class:`.VectorSpaceBasis` (or
          :class:`.MixedVectorSpaceBasis`) spanning the null space of
@@ -224,7 +235,7 @@ def _la_solve(A, x, b, **kwargs):
 
         _la_solve(A, x, b, solver_parameters=parameters_dict)."""
 
-    bcs, solver_parameters, nullspace, nullspace_T, near_nullspace, \
+    P, bcs, solver_parameters, nullspace, nullspace_T, near_nullspace, \
         options_prefix = _extract_linear_solver_args(A, x, b, **kwargs)
 
     # Check whether solution is valid
@@ -234,7 +245,7 @@ def _la_solve(A, x, b, **kwargs):
     if bcs is not None:
         raise RuntimeError("It is no longer possible to apply or change boundary conditions after assembling the matrix `A`; pass any necessary boundary conditions to `assemble` when assembling `A`.")
 
-    solver = ls.LinearSolver(A, solver_parameters=solver_parameters,
+    solver = ls.LinearSolver(A=A, P=P, solver_parameters=solver_parameters,
                              nullspace=nullspace,
                              transpose_nullspace=nullspace_T,
                              near_nullspace=near_nullspace,
@@ -257,7 +268,7 @@ def _la_solve(A, x, b, **kwargs):
 
 
 def _extract_linear_solver_args(*args, **kwargs):
-    valid_kwargs = ["bcs", "solver_parameters", "nullspace",
+    valid_kwargs = ["P", "bcs", "solver_parameters", "nullspace",
                     "transpose_nullspace", "near_nullspace", "options_prefix"]
     if len(args) != 3:
         raise RuntimeError("Missing required arguments, expecting solve(A, x, b, **kwargs)")
@@ -267,14 +278,15 @@ def _extract_linear_solver_args(*args, **kwargs):
             raise RuntimeError("Illegal keyword argument '%s'; valid keywords are %s" %
                                (kwarg, ", ".join("'%s'" % kw for kw in valid_kwargs)))
 
+    P = kwargs.get("P", None)
     bcs = kwargs.get("bcs", None)
     solver_parameters = kwargs.get("solver_parameters", {})
     nullspace = kwargs.get("nullspace", None)
     nullspace_T = kwargs.get("transpose_nullspace", None)
     near_nullspace = kwargs.get("near_nullspace", None)
     options_prefix = kwargs.get("options_prefix", None)
 
-    return bcs, solver_parameters, nullspace, nullspace_T, near_nullspace, options_prefix
+    return P, bcs, solver_parameters, nullspace, nullspace_T, near_nullspace, options_prefix
 
 
 def _extract_args(*args, **kwargs):

firedrake/supermeshing.py

@@ -431,11 +431,14 @@ def likely(cell_A):
     includes = ["-I%s/include" % d for d in dirs]
     libs = ["-L%s/lib" % d for d in dirs]
     libs = libs + ["-Wl,-rpath,%s/lib" % d for d in dirs] + ["-lpetsc", "-lsupermesh"]
-    lib = load(supermesh_kernel_str, "c", "supermesh_kernel",
-               cppargs=includes,
-               ldargs=libs,
-               argtypes=[ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp],
-               restype=ctypes.c_int)
+    lib = load(
+        supermesh_kernel_str, "c", "supermesh_kernel",
+        cppargs=includes,
+        ldargs=libs,
+        argtypes=[ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp, ctypes.c_voidp],
+        restype=ctypes.c_int,
+        comm=mesh_A._comm
+    )
 
     ammm(V_A, V_B, likely, node_locations_A, node_locations_B, M_SS, ctypes.addressof(lib), mat)
     if orig_value_size == 1:

firedrake/utility_meshes.py

@@ -64,8 +64,8 @@
 ]
 
 
-distribution_parameters_noop = {"partition": False,
-                                "overlap_type": (mesh.DistributedMeshOverlapType.NONE, 0)}
+distribution_parameters_no_overlap = {"partition": True,
+                                      "overlap_type": (mesh.DistributedMeshOverlapType.NONE, 0)}
 reorder_noop = False
 
 
@@ -252,7 +252,7 @@ def PeriodicIntervalMesh(
         )
     m = CircleManifoldMesh(
         ncells,
-        distribution_parameters=distribution_parameters_noop,
+        distribution_parameters=distribution_parameters_no_overlap,
         reorder=reorder_noop,
         comm=comm,
         name=name,
@@ -984,7 +984,7 @@ def PeriodicRectangleMesh(
         0.5,
         quadrilateral=quadrilateral,
         reorder=reorder_noop,
-        distribution_parameters=distribution_parameters_noop,
+        distribution_parameters=distribution_parameters_no_overlap,
         comm=comm,
         name=name,
         distribution_name=distribution_name,
@@ -1905,7 +1905,7 @@ def PeriodicBoxMesh(
     m = mesh.Mesh(
         plex,
         reorder=reorder_noop,
-        distribution_parameters=distribution_parameters_noop,
+        distribution_parameters=distribution_parameters_no_overlap,
         name=name,
         distribution_name=distribution_name,
         permutation_name=permutation_name,
@@ -3085,7 +3085,7 @@ def PartiallyPeriodicRectangleMesh(
         longitudinal_direction="z",
         quadrilateral=quadrilateral,
         reorder=reorder_noop,
-        distribution_parameters=distribution_parameters_noop,
+        distribution_parameters=distribution_parameters_no_overlap,
         diagonal=diagonal,
         comm=comm,
         name=name,

tests/regression/test_adjoint_operators.py

@@ -886,3 +886,16 @@ def test_cofunction_subfunctions_with_adjoint():
     k.block_variable.tlm_value = Constant(1)
     get_working_tape().evaluate_tlm()
     assert taylor_test(J_hat, k, Constant(1.0), dJdm=J.block_variable.tlm_value) > 1.9
+
+
+@pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done
+def test_none_riesz_representation_to_derivative():
+    mesh = UnitIntervalMesh(1)
+    space = FunctionSpace(mesh, "Lagrange", 1)
+    u = Function(space).interpolate(SpatialCoordinate(mesh)[0])
+    J = assemble((u ** 2) * dx)
+    rf = ReducedFunctional(J, Control(u))
+    assert isinstance(rf.derivative(), Function)
+    assert isinstance(rf.derivative(options={"riesz_representation": "H1"}), Function)
+    assert isinstance(rf.derivative(options={"riesz_representation": "L2"}), Function)
+    assert isinstance(rf.derivative(options={"riesz_representation": None}), Cofunction)