Constant in dolfinx and ffcx compilation

dolfinx
1

In dolfin-old when defining a form including a constant, the form was not recompiled when changing its value (I suppose there was a specific signature associated to each constant).

This is not anymore the case in dolfinx. When I redefine a constant, the form is recompiled, which is annoying. Is there a way to avoid this?

I give below a MWE. Of course one can assign the value of the constant, but I cannot do it in my case, because the constant is a function that I call for a parametric analysis, and I do not want to have dolfinx.fem.Constant in the launcher of the parameter analysis.

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_params = {"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,c) 
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)
    
# this does not recompile (but I need the former)
cnst = Constant(mesh,1.) 
for c in np.linspace(0,.2,10):
    cnst.value = c
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)
2

The whole point with using form is that it avoids compilation (and caching) of expressions with the same signature.

I think this is a bug in ufl, as ufl does not renumber constants: ufl/form.py at a7426bd8533f2c819f7f164df9c197e277d058c3 · FEniCS/ufl · GitHub and
ufl/form.py at a7426bd8533f2c819f7f164df9c197e277d058c3 · FEniCS/ufl · GitHub

In old dolfin, Constant was not inheriting from ufl.Constant, but ufl.Coefficient, and thus got renumbered.

I’ll see if I can make a working patch.

1 Like
3

@cmaurini I’ve added a PR to ufl that should fix the issue: Renumber constants for caching when defined in loops by jorgensd · Pull Request #113 · FEniCS/ufl · GitHub

Feel free to test it.

1 Like
4

Great, thanks, the fix does the job!

5

Note however the different behaviour between Constant(mesh.ufl_cell(),c) and Constant(mesh,c)
when using different meshes. Why Constant is defined on the whole mesh and not on the single cell in all the examples? What is the rationale for using Constant(mesh.ufl_cell(),c) instead that Constant(mesh,c)?

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


cache_dir = f"{str(Path.cwd())}/cache"
shutil.rmtree(cache_dir,ignore_errors=True)
jit_params = {"cache_dir": cache_dir, "cffi_verbose": True}

# Does not recompile
for c in np.linspace(0,1,10):
    mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
    cnst = Constant(mesh.ufl_cell(),c) 
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)
    
shutil.rmtree(cache_dir,ignore_errors=True)

# Recompile each time
for c in np.linspace(0,1,10):
    mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
    cnst = Constant(mesh,c) 
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)
6

I guess we haven’t thought about this in the examples. The rationale would be for problems where you use cnst*ufl.dx and not cnst*ufl.dx(domain=mesh), as the former will not compile:

for c in np.linspace(0, 1, 10):
    mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
    cnst = Constant(mesh.ufl_cell(), c)
    f = assemble_scalar(form(cnst * dx), jit_params=jit_params)  # (domain=mesh), jit_params=jit_params))
    print(f)

returns

RuntimeError: Expecting to find a Mesh in the form.

I’ve now updated the PR in UFL to make sure that we now use the correct signature for both the cases above.

7

@cmaurini the issue has now been resolved (if you pull the latest version of ufl main branch)

1 Like
8

Thanks a lot. It works also for different meshes, as I was expecting.

# This does not recompile at each iteration
for c in np.linspace(1,10,10):
    mesh = create_unit_square(MPI.COMM_WORLD, 1, int(c))
    f = assemble_scalar(form(cnst * dx(domain=mesh),jit_params=jit_params))
    print(f)

Have a good day.

1 Like