Closed francesco-ballarin closed 1 year ago
I will fix it today. I was expecting some issues to occur... unfortunately, our tests are not comprehensive enough...
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
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
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).
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
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}")
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]
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
Hi, my RBniCSx nightly run failed tonight after the merge of https://github.com/FEniCS/ffcx/pull/594, with the following error
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.txtIf 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
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.