#594 introduced a regression with certain forms

[francesco-ballarin]

Hi,
my RBniCSx nightly run failed tonight after the merge of #594, with the following error

libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c: In function ‘tabulate_tensor_integral_d0b4b45748c1730ae3264b6fe1db607b5018fa9a’:
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c:813:13: error: lvalue required as decrement operand
  813 | sp_083[3] = --2 + c[1] / c[0];
      |             ^~
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c: In function ‘tabulate_tensor_integral_cc5608c9f28da44e30c58e699b254e2dbcb55cb9’:
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c:1038:13: error: lvalue required as decrement operand
 1038 | sp_083[3] = --2 + c[0] / c[1];
      |             ^~
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c: In function ‘tabulate_tensor_integral_a1775329396972d86b0c0859f417536763fb104a’:
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c:1263:13: error: lvalue required as decrement operand
 1263 | sp_083[3] = --2 + c[0] / c[1];
      |             ^~
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c: In function ‘tabulate_tensor_integral_c4e485b50fad96fee23fe972761cfbebbc93f9a1’:
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c:1488:13: error: lvalue required as decrement operand
 1488 | sp_083[3] = --2 + c[1] / c[0];
      |             ^~

The file libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c is attached (upon renaming to .c.txt, otherwise GitHub complains that the file format is not accepted)
libffcx_forms_a485dd25dfc109abba4e332f7d717e453f8b03b0.c.txt

If the cause of the regression is not immediately evident to you, can you help me reverse engineer which integral is being compiled at the failing lines?
The form I am using in the failing tutorials is very complicated

        # mu is dolfinx.fem.Constant with size 3
        thetas_a: typing.Tuple[  # type: ignore[no-any-unimported]
            typing.Callable[[rbnicsx.backends.SymbolicParameters], ufl.core.expr.Expr], ...
        ] = (
            # subdomains 1 and 7
            lambda mu: - (mu[1] - 2) / mu[0] - (2 * (2 * mu[0] - 2) * (mu[0] - 1)) / (mu[0] * (mu[1] - 2)),
            lambda mu: - mu[0] / (mu[1] - 2),
            lambda mu: - (2 * (mu[0] - 1)) / (mu[1] - 2),
            # subdomains 2 and 8
            lambda mu: 2 - (mu[0] - 1) * (mu[0] - 1) / (mu[1] - 2) - mu[1],
            lambda mu: - 1 / (mu[1] - 2),
            lambda mu: (mu[0] - 1) / (mu[1] - 2),
            # subdomains 3 and 5
            lambda mu: - mu[1] / (mu[0] - 2),
            lambda mu: - (mu[0] - 2) / mu[1] - (2 * (2 * mu[1] - 2) * (mu[1] - 1)) / (mu[1] * (mu[0] - 2)),
            lambda mu: - (2 * (mu[1] - 1)) / (mu[0] - 2),
            # subdomains 4 and 6
            lambda mu: - 1 / (mu[0] - 2),
            lambda mu: 2 - (mu[1] - 1) * (mu[1] - 1) / (mu[0] - 2) - mu[0],
            lambda mu: (mu[1] - 1) / (mu[0] - 2),
            # boundaries 5, 6, 7 and 8
            lambda mu: mu[2]
        )
        operators_a = (
            ufl.inner(u.dx(0), v.dx(0)) * dx(1) + ufl.inner(u.dx(0), v.dx(0)) * dx(7),
            ufl.inner(u.dx(1), v.dx(1)) * dx(1) + ufl.inner(u.dx(1), v.dx(1)) * dx(7),
            (ufl.inner(u.dx(0), v.dx(1)) * dx(1) + ufl.inner(u.dx(1), v.dx(0)) * dx(1)
             - ufl.inner(u.dx(0), v.dx(1)) * dx(7) - ufl.inner(u.dx(1), v.dx(0)) * dx(7)),
            # subdomains 2 and 8
            ufl.inner(u.dx(0), v.dx(0)) * dx(2) + ufl.inner(u.dx(0), v.dx(0)) * dx(8),
            ufl.inner(u.dx(1), v.dx(1)) * dx(2) + ufl.inner(u.dx(1), v.dx(1)) * dx(8),
            (ufl.inner(u.dx(0), v.dx(1)) * dx(2) + ufl.inner(u.dx(1), v.dx(0)) * dx(2)
             - ufl.inner(u.dx(0), v.dx(1)) * dx(8) - ufl.inner(u.dx(1), v.dx(0)) * dx(8)),
            # subdomains 3 and 5
            ufl.inner(u.dx(0), v.dx(0)) * dx(3) + ufl.inner(u.dx(0), v.dx(0)) * dx(5),
            ufl.inner(u.dx(1), v.dx(1)) * dx(3) + ufl.inner(u.dx(1), v.dx(1)) * dx(5),
            (ufl.inner(u.dx(0), v.dx(1)) * dx(3) + ufl.inner(u.dx(1), v.dx(0)) * dx(3)
             - ufl.inner(u.dx(0), v.dx(1)) * dx(5) - ufl.inner(u.dx(1), v.dx(0)) * dx(5)),
            # subdomains 4 and 6
            ufl.inner(u.dx(0), v.dx(0)) * dx(4) + ufl.inner(u.dx(0), v.dx(0)) * dx(6),
            ufl.inner(u.dx(1), v.dx(1)) * dx(4) + ufl.inner(u.dx(1), v.dx(1)) * dx(6),
            (ufl.inner(u.dx(0), v.dx(1)) * dx(4) + ufl.inner(u.dx(1), v.dx(0)) * dx(4)
             - ufl.inner(u.dx(0), v.dx(1)) * dx(6) - ufl.inner(u.dx(1), v.dx(0)) * dx(6)),
            # boundaries 5, 6, 7 and 8
            ufl.inner(u, v) * ds(5) + ufl.inner(u, v) * ds(6) + ufl.inner(u, v) * ds(7) + ufl.inner(u, v) * ds(8)
        )
        a = sum([theta_a(mu_symb) * operator_a for (theta_a, operator_a) in zip(thetas_a, operators_a)])
        dolfinx.fem.form(a)  # FAILS

I have attempted to simplify it to a minimally reproducible example, but any simpler example I could come up with (mixture of dx and ds, mixture of partial deriatives, coefficients with addends with a mix of positive/negative signs) did NOT show the issue.

I am quite confident the issue was introduced in #594 because CI ran ok yesterday night.

[chrisrichardson]

I will fix it today. I was expecting some issues to occur... unfortunately, our tests are not comprehensive enough...

[chrisrichardson]

This seems to be a minimal failing example:

mu = Constant(mesh, shape=(3,))
theta = - (mu[1] - 2) / mu[0] - (2 * (2 * mu[0] - 2) * (mu[0] - 1)) / (mu[0] * (mu[1] - 2))
a = theta * inner(grad(u), grad(v)) * dx

[jorgensd]

Also have issues with this PR with DOLFINX_MPC. Workflow ran yesterday, now yielding: https://github.com/jorgensd/dolfinx_mpc/actions/runs/6082097061/job/16499146654

[jorgensd]

This code also does not run the with the newest version of FFCx https://github.com/jorgensd/dolfinx-tutorial/blob/release/chapter2/hyperelasticity.py (Newton solver fails to converge, so I guess similar issue to what happens in MPC).

[chrisrichardson]

Thanks @jorgensd - I was expecting some failures, really. It is showing up inadequacies in the ffcx CI, which I will plug. I'll take a look at your failing examples.
I have now fixed the one for @francesco-ballarin in the branch https://github.com/FEniCS/ffcx/tree/chris/fix-issue-594

[jorgensd]

Code that runs on 0.6.0-r1 that works but fails with ffcx main (including the fix for issue 594):

import numpy as np
import ufl
from dolfinx import fem, log, mesh
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc

L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))


# -

# We create two python functions for determining the facets to apply boundary conditions to

# +
def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], L)


fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# -

# Next, we create a  marker based on these two functions

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# We then create a function for supplying the boundary condition on the left side, which is fixed.

u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)

# To apply the boundary condition, we identity the dofs located on the facets marked by the `MeshTag`.

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

# Next, we define the body force on the reference configuration (`B`), and nominal (first Piola-Kirchhoff) traction (`T`).

B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))

# Define the test and solution functions on the space $V$

v = ufl.TestFunction(V)
u = fem.Function(V)

# Define kinematic quantities used in the problem

# +
# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# -

# Define the elasticity model via a stored strain energy density function $\psi$, and create the expression for the first Piola-Kirchhoff stress:

# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(domain, PETSc.ScalarType(E / (2 * (1 + nu))))
lmbda = fem.Constant(domain, PETSc.ScalarType(E * nu / ((1 + nu) * (1 - 2 * nu))))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

# ```{admonition} Comparison to linear elasticity
# To illustrate the difference between linear and hyperelasticity, the following lines can be uncommented to solve the linear elasticity problem.
# ```

# +
# P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
# -

# Define the variational form with traction integral over all facets with value 2. We set the quadrature degree for the integrals to 4.

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

# As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.

problem = NonlinearProblem(F, u, bcs)

# and then create and customize the Newton solver

# +
solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"


# Finally, we solve the problem over several time steps, updating the y-component of the traction

log.set_log_level(log.LogLevel.INFO)
tval0 = -1.5
for n in range(1, 10):
    T.value[2] = n * tval0
    num_its, converged = solver.solve(u)
    assert (converged)
    u.x.scatter_forward()
    print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")

[chrisrichardson]

How does it fail? I am getting:

Time step 1, Number of iterations 8, Load [ 0.   0.  -1.5]
Time step 2, Number of iterations 9, Load [ 0.  0. -3.]
Time step 3, Number of iterations 10, Load [ 0.   0.  -4.5]
Time step 4, Number of iterations 9, Load [ 0.  0. -6.]
Time step 5, Number of iterations 8, Load [ 0.   0.  -7.5]
Time step 6, Number of iterations 7, Load [ 0.  0. -9.]
Time step 7, Number of iterations 6, Load [  0.    0.  -10.5]
Time step 8, Number of iterations 6, Load [  0.   0. -12.]
Time step 9, Number of iterations 6, Load [  0.    0.  -13.5]

[jorgensd]

How does it fail? I am getting:

Time step 1, Number of iterations 8, Load [ 0.   0.  -1.5]

Time step 2, Number of iterations 9, Load [ 0.  0. -3.]

Time step 3, Number of iterations 10, Load [ 0.   0.  -4.5]

Time step 4, Number of iterations 9, Load [ 0.  0. -6.]

Time step 5, Number of iterations 8, Load [ 0.   0.  -7.5]

Time step 6, Number of iterations 7, Load [ 0.  0. -9.]

Time step 7, Number of iterations 6, Load [  0.    0.  -10.5]

Time step 8, Number of iterations 6, Load [  0.   0. -12.]

Time step 9, Number of iterations 6, Load [  0.    0.  -13.5]

I Get that the newton solver never reaches convergence. This happens on my CI (and locally) https://github.com/jorgensd/dolfinx-tutorial/actions/runs/6082502881/job/16500423176