Fix signature for forms with identitcal constants and identical meshes (up to ufl_id)

[jorgensd]

Reported at: https://fenicsproject.discourse.group/t/constant-in-dolfinx-and-ffcx-compilation/8968?u=dokken

Minimal example failing on main:

import ufl

mesh = ufl.Mesh(ufl.VectorElement("Lagrange", ufl.triangle, 1))
c = ufl.Constant(mesh)
d = ufl.Constant(mesh)
a = c * ufl.dx
b = d * ufl.dx
assert(a.signature() == b.signature())

Constants now has the same behavior as coefficients:

V = ufl.FunctionSpace(mesh, ufl.FiniteElement("CG", mesh.ufl_cell(), 1))
u = ufl.Coefficient(V)
v = ufl.TestFunction(V)
w = ufl.Coefficient(V)

a = ufl.inner(u, v) * ufl.dx
d = ufl.inner(w, v) * ufl.dx
assert(a.signature() == d.signature())

This PR now also supports same signature for different meshes (with same coordinate element), i.e.

mesh = ufl.Mesh(ufl.VectorElement("Lagrange", ufl.triangle, 1))
mesh2 = ufl.Mesh(ufl.VectorElement("Lagrange", ufl.triangle, 1))
c = ufl.Constant(mesh)
d = ufl.Constant(mesh2)
a = c * ufl.dx
b = d * ufl.dx
assert(a.signature() == b.signature())

[jorgensd]

@dham This change should not affect firedrake, as it is using ufl.Coefficient for constants (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/constant.py#L26) which is the same as legacy dolfin does.

[dham]

@dham This change should not affect firedrake, as it is using ufl.Coefficient for constants (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/constant.py#L26) which is the same as legacy dolfin does.

I'm sure that's right. I've just set https://github.com/firedrakeproject/firedrake/pull/2521 running to be sure. If that passes then I'm happy.

[twmr]

There still seems to be problems in builds with a complex PetscScalar.

I've tried running the demo from the linked discourse issue. What I observe is that when I pass a float to fem.Constant it fails to compile the form, whereas if I cast the number to a complex, it works.

However, for this part of the MWE

# this does not recompile (but I need the former)
cnst = Constant(mesh,1.) 
for c in np.linspace(0,.2,10):
    cnst.value = complex(c, 0)
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)

It doesnt even work when I pass a complex to cnst.value.

Traceback (most recent call last):
  File "/home/thomas/gitrepos/fenics_laser_states_2022/demos/fenics_qa_ffc_demo.py", line 24, in <module>
    cnst.value = complex(c, 0)
  File "/home/thomas/miniconda/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/function.py", line 63, in value
    np.copyto(self._cpp_object.value, np.asarray(v))
  File "<__array_function__ internals>", line 180, in copyto
TypeError: Cannot cast scalar from dtype('complex128') to dtype('float64') according to the rule 'same_kind'

[jorgensd]

There still seems to be problems in builds with a complex PetscScalar.

I've tried running the demo from the linked discourse issue. What I observe is that when I pass a float to fem.Constant it fails to compile the form, whereas if I cast the number to a complex, it works.

However, for this part of the MWE

# this does not recompile (but I need the former)
cnst = Constant(mesh,1.) 
for c in np.linspace(0,.2,10):
    cnst.value = complex(c, 0)
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)

It doesnt even work when I pass a complex to cnst.value.

Traceback (most recent call last):
  File "/home/thomas/gitrepos/fenics_laser_states_2022/demos/fenics_qa_ffc_demo.py", line 24, in <module>
    cnst.value = complex(c, 0)
  File "/home/thomas/miniconda/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/function.py", line 63, in value
    np.copyto(self._cpp_object.value, np.asarray(v))
  File "<__array_function__ internals>", line 180, in copyto
TypeError: Cannot cast scalar from dtype('complex128') to dtype('float64') according to the rule 'same_kind'

As seen in the source code of DOLFINx: https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/fem/function.py#L32-L57 the constant is created wrt. to the dtype of the original input. I.e. for complex builds you should send in the right dtype at compilation time, see:

from petsc4py import PETSc
from dolfinx.fem import (Constant, assemble_scalar, form)
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
from ufl import dx
import numpy as np
from pathlib import Path
import shutil

mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)

cache_dir = f"{str(Path.cwd())}/cache"
shutil.rmtree(cache_dir, ignore_errors=True)
jit_options = {"cache_dir": cache_dir, "cffi_verbose": True}
# this recompiles the form at each iterations (different behavior wrt dolfin-old)
for c in np.linspace(0, 1, 10):
    cnst = Constant(mesh, PETSc.ScalarType(c))
    f = assemble_scalar(form(cnst * dx(domain=mesh), jit_options=jit_options))
    print(f)

# this does not recompile (but I need the former)
cnst = Constant(mesh, PETSc.ScalarType(1.))
for c in np.linspace(0, .2, 10):
    cnst.value = c
    f = assemble_scalar(form(cnst * dx(domain=mesh), jit_options=jit_options))
    print(f)

for an example compatible with dolfinx/dolfinx:nightly built the 20th of October. The only difference for the 0.5.1 release is the change of jit_params to jit_options.

[twmr]

Thx @jorgensd, this works!