Fix bug in coefficient packing

[garth-wells]

Packing got all coefficients for an integral (domain) type, rather than just coefficients that are present in an integral over a subdomain (i.e., an integral with an ID).

[jorgensd]

Is there a test we could add to capture this subtle change?

[jorgensd]

For instance, main yields:

 end-start=4.45428e-01
[[  1.   1.   1. ... 300. 300. 300.]]

while your branch yields

end-start=7.89978e-02
[[0. 0. 0. ... 0. 0. 0.]]

running

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
M = 5
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, M,M,M)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
exterior_facet_indices = dolfinx.mesh.exterior_facet_indices(mesh.topology)
facet_markers= np.arange(len(exterior_facet_indices), dtype=np.int32)
ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, exterior_facet_indices, facet_markers)

N = len(exterior_facet_indices)
us = [dolfinx.fem.Function(V) for _ in range(N)]
for i,u in enumerate(us):
    u.x.array[:] = i + 1
    u.name=f"u_{i}"

ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)
volume_form = us[0] * ufl.dx
surface_form = sum(us[i] * ds(i) for i in range(N))

combined_form = volume_form + surface_form
compiled_form = dolfinx.fem.form(combined_form)

import time
start = time.perf_counter()
coeffs = dolfinx.fem.assemble.pack_coefficients(compiled_form._cpp_object)
end =time.perf_counter()
print(f"{dolfinx.common.git_commit_hash=} {end-start=:.5e}")
print(coeffs[(dolfinx.fem.IntegralType.exterior_facet, 1)])

[garth-wells]

For instance, main yields:

 end-start=4.45428e-01
[[  1.   1.   1. ... 300. 300. 300.]]

while your branch yields

end-start=7.89978e-02
[[0. 0. 0. ... 0. 0. 0.]]

Which is correct?

[jorgensd]

For instance, main yields:

 end-start=4.45428e-01
[[  1.   1.   1. ... 300. 300. 300.]]

while your branch yields

end-start=7.89978e-02
[[0. 0. 0. ... 0. 0. 0.]]

Which is correct?

Yes, you can verify this by

For instance, main yields:

 end-start=4.45428e-01
[[  1.   1.   1. ... 300. 300. 300.]]

while your branch yields

end-start=7.89978e-02
[[0. 0. 0. ... 0. 0. 0.]]

Which is correct?

Yes. You can verify the correctness by for instance checking for a single integral id, here (1):

coeff = coeffs[(dolfinx.fem.IntegralType.exterior_facet, 1)]
print(np.flatnonzero(coeff.flatten()), coeff[0,np.flatnonzero(coeff.flatten())])

which yields:

[[0. 0. 0. ... 0. 0. 0.]]
[20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39] [2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]

which makes sense as:

1. The integral with subdomain_id=1 contains a single coefficient (and a single cell).
2. There are 20 DOFS in a 3rd order Lagrange space, meaning that the first coefficient (related to integral id=0) is not packed, while the next one is packed. All others are 0 as they are not part of u[1]*ds(1).

[garth-wells]

#3646 includes this change, and triggers an error if coefficients that are not required by an integral are packed.

It triggers because in #3646 the Form will not have any entity data to pack over if the coefficient is not required by the form - it worked until now because Form would compute data on-demand (repeatedly) that is not required.