Slow Newton convergence using dolfinx_materials

[JakobDegn]

Hello

This is my first post here so please let me know if anything is not following the current guidelines. I have tried follow the post (Read before posting: How do I get my question answered?) the best I could.

Anyways, I am having trouble with slow convergence of a mechanical problem defined with the library dolfinx_materials. The problem concerns a simple linear elastic bar in 2 dimensions which is fixed in one end and has a prescribed displacement in the other end. For simplicity, the mesh consists of only 2 linear quadrilateral elements.

The behavior of the bar is defined in MFront using the volumetric and shear components, i.e.

\boldsymbol{\sigma} = \kappa \operatorname{tr} [\boldsymbol{\varepsilon}] \boldsymbol{I} + 2 \mu \operatorname{dev} [\boldsymbol{\varepsilon}]

\mathbb{C} = \kappa \; \boldsymbol{I} \otimes \boldsymbol{I} + 2 \mu \mathbb{P}_\text{dev}

Using the following scripts, the problem converges in 8 iterations. Comparatively, if I define the residual and Jacobian using only UFL (i.e. without MFront and dolfinx_materials), the problem converges in 1 iteration as expected.

Scripts

The FEniCSx script:

############################################
# PREAMBLE
############################################

# Import python libraries
import numpy as np
import os
from petsc4py import PETSc
from mpi4py import MPI

# Import fenicsx libraries
from dolfinx import fem, mesh
from dolfinx.fem import petsc
from basix.ufl import element
import ufl

# Import dolfinx_materials libraries
from dolfinx_materials.quadrature_map import QuadratureMap
from dolfinx_materials.material.mfront import MFrontMaterial
from dolfinx_materials.solvers import SNESNonlinearMaterialProblem
from dolfinx_materials.utils import symmetric_tensor_to_vector

############################################
# MESH 
############################################

# Parallel communicator
comm = MPI.COMM_WORLD

# Create the mesh
msh = mesh.create_rectangle(comm=comm, points=[[0.0, 0.0], [25.0, 1.0]], n=[2, 1], cell_type=mesh.CellType.quadrilateral)

# Dimensionality
tdim = msh.topology.dim
fdim = tdim - 1

# Get boundary edges
edge_right = mesh.locate_entities_boundary(mesh=msh, dim=fdim, marker=lambda x: np.isclose(x[0], 25.0))
edge_left = mesh.locate_entities_boundary(mesh=msh, dim=fdim, marker=lambda x: np.isclose(x[0], 0.0))

############################################
# FUNCTIONS AND SPACES
############################################

# Define element
ele = element(family="Lagrange", cell=msh.basix_cell(), degree=1, shape=(tdim,))

# Create function space
V = fem.functionspace(mesh=msh, element=ele)

# Define current state function
u = fem.Function(V, name="displacement")

# Define test/trial functions
v = ufl.TestFunction(V)
du = ufl.TrialFunction(V)

############################################
# BOUNDARY CONDITIONS
############################################

# Get boundary displacement dofs
dofs_ux_right = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=edge_right)
dofs_ux_left = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=edge_left)
dofs_uy_left = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=edge_left)

# Set displacement boundary condition values
u_load = fem.Constant(domain=msh, c=PETSc.ScalarType(1.0))
u_fixed = fem.Constant(domain=msh, c=PETSc.ScalarType(0.0))

# Define Dirichlet BCs for displacement
bc_ux_right = fem.dirichletbc(value=u_load, dofs=dofs_ux_right, V=V.sub(0))
bc_ux_left = fem.dirichletbc(value=u_fixed, dofs=dofs_ux_left, V=V.sub(0))
bc_uy_left = fem.dirichletbc(value=u_fixed, dofs=dofs_uy_left, V=V.sub(1))

# Collect boundary conditions
bcs = [bc_ux_right, bc_ux_left, bc_uy_left]

############################################
# MATERIAL DEFINITION
############################################

# Define MFront material
material = MFrontMaterial(
    path=os.path.join(os.getcwd(), "mfront/src/libBehaviour.so"),
    name="DisplacementDeviatoric",
    hypothesis="plane_strain",
    material_properties={"YoungModulus": 1.0, "PoissonRatio": 0.0},
)

# Define quadrature map
qmap = QuadratureMap(mesh=msh, deg=2, material=material)

############################################
# VARIATIONAL FORMULATION
############################################

# Displacement functions
def eps(u):
    '''Strain tensor as vector'''
    return symmetric_tensor_to_vector(ufl.sym(ufl.grad(u)))

# Declare strain measure
qmap.register_gradient("Strain", eps(u))

# Extract stress from material
sig = qmap.fluxes["Stress"]

# Displacement subproblem
res = ufl.dot(sig, eps(v)) * qmap.dx
Jac = qmap.derivative(res, u, du)

############################################
# SOLVER
############################################

# Define variational problem
problem = SNESNonlinearMaterialProblem(qmap=qmap, F=res, J=Jac, u=u, bcs=bcs)

# Create and configure solver
solver = PETSc.SNES().create(comm)
solver.setType("newtonls")
solver.setTolerances(rtol=0.0, atol=1e-6, stol=1e-8, max_it=200)
solver.getKSP().setType("preonly")
solver.getKSP().setTolerances(rtol=1.0e-9)
solver.getKSP().getPC().setType("lu")
solver.getKSP().getPC().setFactorSolverType('mumps')

# Disable line search
PETSc.Options().setValue("snes_linesearch_type", "basic")
solver.setFromOptions()

############################################
# TEST SOLVER
############################################

# Set initial guess
petsc.set_bc(u.x.petsc_vec, bcs)
u.x.scatter_forward()

# Solve
problem.solve(solver, print_solution=True)

The MFront behavior:

@DSL DefaultDSL; 
@Behaviour DisplacementDeviatoric;

// ╔════════════╗
// ║ Hypotheses ║
// ╚════════════╝

@ModellingHypothesis PlaneStrain;

// ╔═════════════════════╗
// ║ Material properties ║
// ╚═════════════════════╝

@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");

// ╔═════════════════╗
// ║ Local Variables ║
// ╚═════════════════╝

@LocalVariable stress λ; // first Lamé parameter
@LocalVariable stress μ; // second Lamé parameter (shear modulus)
@LocalVariable stress κ; // bulk modulus

// ╔══════════════════════╗
// ║ Init Local Variables ║
// ╚══════════════════════╝

@InitLocalVariables {
  λ = computeLambda(young, nu);
  μ = computeMu(young, nu);
  κ = λ + (2. / 3.) * μ;
}

// ╔════════════╗
// ║ Integrator ║
// ╚════════════╝

@ProvidesSymmetricTangentOperator;
@Integrator {
  // current total strain tensor
  const auto ε = eval(εᵗᵒ + Δεᵗᵒ);
  // stress computation
  σ = κ * trace(ε) * I₂ + 2 * μ * deviator(ε);
  // tangent operator
  if (computeTangentOperator_) {
    // suppress unused variable warning 
    static_cast<void>(smt);
    // compute tangent operator
    ∂σ∕∂Δεᵗᵒ = κ * Stensor4::IxI() + 2 * μ * Stensor4::K();
  } // end if computeTangentOperator_
} // end of @Integrator

The pure UFL implementation is practically identical to the dolfinx_materials implementation except for the residual and Jacobian which are defined like this:

# Material parameters
lmbda = E * nu / (1 - nu**2)
mu = E / (2 * (1 + nu))
kappa = lmbda + 2/tdim*mu

# Stress and strain tensors
eps = ufl.variable(ufl.sym(ufl.grad(u)))
sig = kappa * ufl.tr(eps) * ufl.Identity(tdim) + 2 * mu * ufl.dev(eps)

# Problem definition
res = ufl.inner(sig, ufl.grad(v)) * ufl.dx
Jac = ufl.derivative(res, u, du)

Additional information

Convergence history

This is the convergence history of residual norms:

dolfinx_materials implementation:

iter = 0, ||r|| = 5.656854e-02
iter = 1, ||r|| = 1.118034e-02
iter = 2, ||r|| = 2.795085e-03
iter = 3, ||r|| = 6.987712e-04
iter = 4, ||r|| = 1.746928e-04
iter = 5, ||r|| = 4.367320e-05
iter = 6, ||r|| = 1.091830e-05
iter = 7, ||r|| = 2.729575e-06
iter = 8, ||r|| = 6.823938e-07

pure UFL implementation:

iter = 0, ||r|| = 5.656854e-02
iter = 1, ||r|| = 1.182172e-16

Which suggests that at least the residual is computed correctly as the initial residuals are identical.

Jacobian matrix

The Jacobian matrix for the two implementations:

dolfinx_materials implementation:

[[ 1.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    1.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    1.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    1.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    5.61  1.96 -5.53 -1.96  0.    0.32  0.   -0.66]
 [ 0.    0.    0.    0.    0.01  8.37  0.01 -8.32  0.    2.12  0.   -2.15]
 [ 0.    0.    0.    0.   -5.53 -1.96  5.61  1.96  0.   -0.32  0.    0.66]
 [ 0.    0.    0.    0.    0.01 -8.32  0.01  8.37  0.   -2.03  0.    2.01]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    1.    0.    0.    0.  ]
 [ 0.    0.    0.    0.   -0.17  2.01  0.16 -2.03  0.    4.13  0.   -4.1 ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    1.    0.  ]
 [ 0.    0.    0.    0.   -0.17 -2.15  0.16  2.12  0.   -4.22  0.    4.24]]

pure UFL implementation:

[[ 1.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    1.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    1.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    1.    0.    0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    4.22 -0.   -4.14 -0.    0.   -0.12  0.   -0.12]
 [ 0.    0.    0.    0.   -0.    8.36  0.   -8.32  0.    2.07  0.   -2.09]
 [ 0.    0.    0.    0.   -4.14  0.    4.22 -0.    0.    0.12  0.    0.12]
 [ 0.    0.    0.    0.   -0.   -8.32 -0.    8.36  0.   -2.09  0.    2.07]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    1.    0.    0.    0.  ]
 [ 0.    0.    0.    0.   -0.12  2.07  0.12 -2.09  0.    4.18  0.   -4.16]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    1.    0.  ]
 [ 0.    0.    0.    0.   -0.12 -2.09  0.12  2.07  0.   -4.16  0.    4.18]]

The main thing I notice is that the Jacobian of the dolfinx_materials implementation is not symmetric. I imagine there is a (stupid) error somewhere in my script that I just can’t find. Otherwise I think my problem might have something to do with the way dolfinx_materials define the Jacobian?

Installation and platform

I have installed FEniCSx from conda-forge using micromamba. The FEniCSx version is 0.8.0 and I have cloned the main branch of dolfinx_materials I am working in VS Code with a DevContainer setup, but honestly I do not know a lot about this, so anyone thinks this could influence the convergence issue i am having, I will be happy to share the configuration files.

I am grateful for any advice

[JakobDegn]

If anybody comes across the same issue, i found a solution to my problem.

It seems like there happens something internally to the MFront function Stensor4::IxI(). Defining it before using it apparently solves this problem, i.e. the @Integrator block in my MFront file becomes

// ╔════════════╗
// ║ Integrator ║
// ╚════════════╝

@ProvidesSymmetricTangentOperator;
@Integrator {
  // current total strain tensor
  const auto ε = eval(εᵗᵒ + Δεᵗᵒ);
  // stress computation
  σ = κ * trace(ε) * I₂ + 2 * μ * deviator(ε);
  // Fourth-order tensor: I₂⊗I₂
  constexpr Stensor4 IxI = Stensor4::IxI();
  // tangent operator
  if (computeTangentOperator_) {
    // suppress unused variable warning 
    static_cast<void>(smt);
    // compute tangent operator
    ∂σ∕∂Δεᵗᵒ = κ * IxI + 2 * μ * Stensor4::K();
  } // end if computeTangentOperator_
} // end of @Integrator

With this definition, my problem now converges in 1 iteration.