#594 introduced a regression with certain forms
Closed this issue · 8 comments
francesco-ballarin commented
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 commented
I will fix it today. I was expecting some issues to occur... unfortunately, our tests are not comprehensive enough...
chrisrichardson commented
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 commented
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 commented
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 commented
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 commented
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 commented
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 commented
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