Plasticity of Mohr-Coulomb with apex-smoothing (JAX)#
This tutorial aims to demonstrate how modern automatic algorithmic
differentiation (AD) techniques may be used to define a complex constitutive
model demanding a lot of by-hand differentiation. In particular, we implement
the non-associative plasticity model of Mohr-Coulomb with apex-smoothing applied
to a slope stability problem for soil. We use the
JAX package to define constitutive
relations including the differentiation of certain terms and
FEMExternalOperator class to incorporate this model into a weak formulation
within UFL.
The tutorial is based on the limit analysis within semi-definite programming framework, where the plasticity model was replaced by the MFront/TFEL implementation of the Mohr-Coulomb elastoplastic model with apex smoothing.
Problem formulation#
We solve a slope stability problem of a soil domain \(\Omega\) represented by a rectangle \([0; L] \times [0; W]\) with homogeneous Dirichlet boundary conditions for the displacement field \(\boldsymbol{u} = \boldsymbol{0}\) on the right side \(x = L\) and the bottom one \(z = 0\). The loading consists of a gravitational body force \(\boldsymbol{q}=[0, -\gamma]^T\) with \(\gamma\) being the soil self-weight. The solution of the problem is to find the collapse load \(q_\text{lim}\), for which we know an analytical solution in the case of the standard Mohr-Coulomb model without smoothing under plane strain assumption for associative plastic law [Chen and Liu, 1990]. Here we follow the same Mandel-Voigt notation as in the von Mises plasticity tutorial.
If \(V\) is a functional space of admissible displacement fields, then we can write out a weak formulation of the problem:
Find \(\boldsymbol{u} \in V\) such that
\[
F(\boldsymbol{u}; \boldsymbol{v}) = \int\limits_\Omega
\boldsymbol{\sigma}(\boldsymbol{u}) \cdot
\boldsymbol{\varepsilon}(\boldsymbol{v}) \, \mathrm{d}\boldsymbol{x} -
\int\limits_\Omega \boldsymbol{q} \cdot \boldsymbol{v} \, \mathrm{d}\boldsymbol{x} = \boldsymbol{0}, \quad
\forall \boldsymbol{v} \in V,
\]
where \(\boldsymbol{\sigma}\) is an external operator representing the stress tensor.
Note
Although the tutorial shows the implementation of the Mohr-Coulomb model, it is quite general to be adapted to a wide rage of plasticity models that may be defined through a yield surface and a plastic potential.
Implementation#
Preamble#
from functools import partial
from mpi4py import MPI
from petsc4py import PETSc
import jax
import jax.lax
import jax.numpy as jnp
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
from mpltools import annotation # for slope markers
from utilities import find_cell_by_point
import basix
import ufl
from dolfinx import default_scalar_type, fem, mesh
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators,
)
from dolfinx_external_operator.petsc import assemble_residual_with_callback
jax.config.update("jax_enable_x64", True)
Here we define geometrical and material parameters of the problem as well as some useful constants.
E = 6778 # [MPa] Young modulus
nu = 0.25 # [-] Poisson ratio
c = 3.45 # [MPa] cohesion
phi = 30 * np.pi / 180 # [rad] friction angle
psi = 30 * np.pi / 180 # [rad] dilatancy angle
theta_T = 26 * np.pi / 180 # [rad] transition angle as defined by Abbo and Sloan
a = 0.26 * c / np.tan(phi) # [MPa] tension cuff-off parameter
L, H = (1.2, 1.0)
Nx, Ny = (25, 25)
gamma = 1.0
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([L, H])], [Nx, Ny])
k_u = 2
gdim = domain.topology.dim
V = fem.functionspace(domain, ("Lagrange", k_u, (gdim,)))
# Boundary conditions
def on_right(x):
return np.isclose(x[0], L)
def on_bottom(x):
return np.isclose(x[1], 0.0)
bottom_dofs = fem.locate_dofs_geometrical(V, on_bottom)
right_dofs = fem.locate_dofs_geometrical(V, on_right)
bcs = [
fem.dirichletbc(np.array([0.0, 0.0], dtype=PETSc.ScalarType), bottom_dofs, V),
fem.dirichletbc(np.array([0.0, 0.0], dtype=PETSc.ScalarType), right_dofs, V),
]
def epsilon(v):
grad_v = ufl.grad(v)
return ufl.as_vector(
[
grad_v[0, 0],
grad_v[1, 1],
0.0,
np.sqrt(2.0) * 0.5 * (grad_v[0, 1] + grad_v[1, 0]),
]
)
k_stress = 2 * (k_u - 1)
dx = ufl.Measure(
"dx",
domain=domain,
metadata={"quadrature_degree": k_stress, "quadrature_scheme": "default"},
)
stress_dim = 2 * gdim
S_element = basix.ufl.quadrature_element(domain.topology.cell_name(), degree=k_stress, value_shape=(stress_dim,))
S = fem.functionspace(domain, S_element)
Du = fem.Function(V, name="Du")
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="du")
v = ufl.TestFunction(V)
sigma = FEMExternalOperator(epsilon(Du), function_space=S)
sigma_n = fem.Function(S, name="sigma_n")
Defining plasticity model and external operator#
The constitutive model of the soil is described by a non-associative plasticity law without hardening that is defined by the Mohr-Coulomb yield surface \(f\) and the plastic potential \(g\). Both quantities may be expressed through the following function \(h\)
\[\begin{align*}
& h(\boldsymbol{\sigma}, \alpha) =
\frac{I_1(\boldsymbol{\sigma})}{3}\sin\alpha +
\sqrt{J_2(\boldsymbol{\sigma}) K^2(\alpha) + a^2(\alpha)\sin^2\alpha} -
c\cos\alpha, \\
& f(\boldsymbol{\sigma}) = h(\boldsymbol{\sigma}, \phi), \\
& g(\boldsymbol{\sigma}) = h(\boldsymbol{\sigma}, \psi),
\end{align*}\]
where \(\phi\) and \(\psi\) are friction and dilatancy angles, \(c\) is a cohesion, \(I_1(\boldsymbol{\sigma}) = \mathrm{tr} \boldsymbol{\sigma}\) is the first invariant of the stress tensor and \(J_2(\boldsymbol{\sigma}) = \frac{1}{2}\boldsymbol{s} \cdot \boldsymbol{s}\) is the second invariant of the deviatoric part of the stress tensor. The expression of the coefficient \(K(\alpha)\) may be found in the MFront/TFEL implementation of this plastic model.
During the plastic loading the stress-strain state of the solid must satisfy the following system of nonlinear equations
(5)#\[\begin{split} \begin{cases}
\boldsymbol{r}_{g}(\boldsymbol{\sigma}_{n+1}, \Delta\lambda) =
\boldsymbol{\sigma}_{n+1} - \boldsymbol{\sigma}_n -
\boldsymbol{C} \cdot (\Delta\boldsymbol{\varepsilon} - \Delta\lambda
\frac{\mathrm{d} g}{\mathrm{d}\boldsymbol{\sigma}}(\boldsymbol{\sigma_{n+1}})) =
\boldsymbol{0}, \\
r_f(\boldsymbol{\sigma}_{n+1}) = f(\boldsymbol{\sigma}_{n+1}) = 0,
\end{cases}\end{split}\]
where \(\Delta\) is associated with increments of a quantity between the next loading step \(n + 1\) and the current loading step \(n\).
By introducing the residual vector \(\boldsymbol{r} = [\boldsymbol{r}_{g}^T, r_f]^T\) and its argument vector \(\boldsymbol{y}_{n+1} = [\boldsymbol{\sigma}_{n+1}^T, \Delta\lambda]^T\), we obtain the following nonlinear constitutive equation:
\[
\boldsymbol{r}(\boldsymbol{y}_{n+1}) = \boldsymbol{0}.
\]
To solve this equation we apply the Newton method and introduce the local Jacobian of the residual vector \(\boldsymbol{j} := \frac{\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{y}}\). Thus we solve the following linear system at each quadrature point for the plastic phase
\[\begin{split}
\begin{cases}
\boldsymbol{j}(\boldsymbol{y}_{n})\boldsymbol{t} = -
\boldsymbol{r}(\boldsymbol{y}_{n}), \\
\boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \boldsymbol{t}.
\end{cases}
\end{split}\]
During the elastic loading, we consider a trivial system of equations
(6)#\[\begin{split}
\begin{cases}
\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_n +
\boldsymbol{C} \cdot \Delta\boldsymbol{\varepsilon}, \\ \Delta\lambda = 0.
\end{cases}
\end{split}\]
The algorithm solving the systems (5)–(6) is called the return-mapping procedure and the solution defines the return-mapping correction of the stress tensor. By implementation of the external operator \(\boldsymbol{\sigma}\) we mean the implementation of this algorithmic procedure.
The automatic differentiation tools of the JAX library are applied to calculate the three distinct derivatives:
\(\frac{\mathrm{d} g}{\mathrm{d}\boldsymbol{\sigma}}\) - derivative of the plastic potential \(g\),
\(j = \frac{\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{y}}\) - derivative of the local residual \(\boldsymbol{r}\),
\(\boldsymbol{C}_\text{tang} = \frac{\mathrm{d}\boldsymbol{\sigma}}{\mathrm{d}\boldsymbol{\varepsilon}}\) - stress tensor derivative or consistent tangent moduli.
Defining yield surface and plastic potential#
First of all, we define supplementary functions that help us to express the
yield surface \(f\) and the plastic potential \(g\). In the following definitions,
we use built-in functions of the JAX package, in particular, the conditional
primitive jax.lax.cond. It is necessary for the correct work of the AD tool
and just-in-time compilation. For more details, please, visit the JAX
documentation.
def J3(s):
return s[2] * (s[0] * s[1] - s[3] * s[3] / 2.0)
def J2(s):
return 0.5 * jnp.vdot(s, s)
def theta(s):
J2_ = J2(s)
arg = -(3.0 * np.sqrt(3.0) * J3(s)) / (2.0 * jnp.sqrt(J2_ * J2_ * J2_))
arg = jnp.clip(arg, -1.0, 1.0)
theta = 1.0 / 3.0 * jnp.arcsin(arg)
return theta
def sign(x):
return jax.lax.cond(x < 0.0, lambda x: -1, lambda x: 1, x)
def coeff1(theta, angle):
return np.cos(theta_T) - (1.0 / np.sqrt(3.0)) * np.sin(angle) * np.sin(theta_T)
def coeff2(theta, angle):
return sign(theta) * np.sin(theta_T) + (1.0 / np.sqrt(3.0)) * np.sin(angle) * np.cos(theta_T)
coeff3 = 18.0 * np.cos(3.0 * theta_T) * np.cos(3.0 * theta_T) * np.cos(3.0 * theta_T)
def C(theta, angle):
return (
-np.cos(3.0 * theta_T) * coeff1(theta, angle) - 3.0 * sign(theta) * np.sin(3.0 * theta_T) * coeff2(theta, angle)
) / coeff3
def B(theta, angle):
return (
sign(theta) * np.sin(6.0 * theta_T) * coeff1(theta, angle) - 6.0 * np.cos(6.0 * theta_T) * coeff2(theta, angle)
) / coeff3
def A(theta, angle):
return (
-(1.0 / np.sqrt(3.0)) * np.sin(angle) * sign(theta) * np.sin(theta_T)
- B(theta, angle) * sign(theta) * np.sin(3 * theta_T)
- C(theta, angle) * np.sin(3.0 * theta_T) * np.sin(3.0 * theta_T)
+ np.cos(theta_T)
)
def K(theta, angle):
def K_false(theta):
return jnp.cos(theta) - (1.0 / np.sqrt(3.0)) * np.sin(angle) * jnp.sin(theta)
def K_true(theta):
return (
A(theta, angle)
+ B(theta, angle) * jnp.sin(3.0 * theta)
+ C(theta, angle) * jnp.sin(3.0 * theta) * jnp.sin(3.0 * theta)
)
return jax.lax.cond(jnp.abs(theta) > theta_T, K_true, K_false, theta)
def a_g(angle):
return a * np.tan(phi) / np.tan(angle)
dev = np.array(
[
[2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0, 0.0],
[-1.0 / 3.0, 2.0 / 3.0, -1.0 / 3.0, 0.0],
[-1.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0, 0.0],
[0.0, 0.0, 0.0, 1.0],
],
dtype=PETSc.ScalarType,
)
tr = np.array([1.0, 1.0, 1.0, 0.0], dtype=PETSc.ScalarType)
def surface(sigma_local, angle):
s = dev @ sigma_local
I1 = tr @ sigma_local
theta_ = theta(s)
return (
(I1 / 3.0 * np.sin(angle))
+ jnp.sqrt(
J2(s) * K(theta_, angle) * K(theta_, angle) + a_g(angle) * a_g(angle) * np.sin(angle) * np.sin(angle)
)
- c * np.cos(angle)
)
By picking up an appropriate angle we define the yield surface \(f\) and the plastic potential \(g\).
def f(sigma_local):
return surface(sigma_local, phi)
def g(sigma_local):
return surface(sigma_local, psi)
dgdsigma = jax.jacfwd(g)
Solving constitutive equations#
In this section, we define the constitutive model by solving the systems
(5)–(6). They must be solved at each Gauss point, so we
apply the Newton method, implement the whole algorithm locally and then
vectorize the final result using jax.vmap.
In the following cell, we define locally the residual \(\boldsymbol{r}\) and
its Jacobian drdy.
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E / (2.0 * (1.0 + nu))
C_elas = np.array(
[
[lmbda + 2 * mu, lmbda, lmbda, 0],
[lmbda, lmbda + 2 * mu, lmbda, 0],
[lmbda, lmbda, lmbda + 2 * mu, 0],
[0, 0, 0, 2 * mu],
],
dtype=PETSc.ScalarType,
)
S_elas = np.linalg.inv(C_elas)
ZERO_VECTOR = np.zeros(stress_dim, dtype=PETSc.ScalarType)
def deps_p(sigma_local, dlambda, deps_local, sigma_n_local):
sigma_elas_local = sigma_n_local + C_elas @ deps_local
yielding = f(sigma_elas_local)
def deps_p_elastic(sigma_local, dlambda):
return ZERO_VECTOR
def deps_p_plastic(sigma_local, dlambda):
return dlambda * dgdsigma(sigma_local)
return jax.lax.cond(yielding <= 0.0, deps_p_elastic, deps_p_plastic, sigma_local, dlambda)
def r_g(sigma_local, dlambda, deps_local, sigma_n_local):
deps_p_local = deps_p(sigma_local, dlambda, deps_local, sigma_n_local)
return sigma_local - sigma_n_local - C_elas @ (deps_local - deps_p_local)
def r_f(sigma_local, dlambda, deps_local, sigma_n_local):
sigma_elas_local = sigma_n_local + C_elas @ deps_local
yielding = f(sigma_elas_local)
def r_f_elastic(sigma_local, dlambda):
return dlambda
def r_f_plastic(sigma_local, dlambda):
return f(sigma_local)
return jax.lax.cond(yielding <= 0.0, r_f_elastic, r_f_plastic, sigma_local, dlambda)
def r(y_local, deps_local, sigma_n_local):
sigma_local = y_local[:stress_dim]
dlambda_local = y_local[-1]
res_g = r_g(sigma_local, dlambda_local, deps_local, sigma_n_local)
res_f = r_f(sigma_local, dlambda_local, deps_local, sigma_n_local)
res = jnp.c_["0,1,-1", res_g, res_f] # concatenates an array and a scalar
return res
drdy = jax.jacfwd(r)
Then we define the function return_mapping that implements the
return-mapping algorithm numerically via the Newton method.
Nitermax, tol = 200, 1e-8
ZERO_SCALAR = np.array([0.0])
def return_mapping(deps_local, sigma_n_local):
"""Performs the return-mapping procedure.
It solves elastoplastic constitutive equations numerically by applying the
Newton method in a single Gauss point. The Newton loop is implement via
`jax.lax.while_loop`.
The function returns `sigma_local` two times to reuse its values after
differentiation, i.e. as once we apply
`jax.jacfwd(return_mapping, has_aux=True)` the ouput function will
have an output of
`(C_tang_local, (sigma_local, niter_total, yielding, norm_res, dlambda))`.
Returns:
sigma_local: The stress at the current Gauss point.
niter_total: The total number of iterations.
yielding: The value of the yield function.
norm_res: The norm of the residuals.
dlambda: The value of the plastic multiplier.
"""
niter = 0
dlambda = ZERO_SCALAR
sigma_local = sigma_n_local
y_local = jnp.concatenate([sigma_local, dlambda])
res = r(y_local, deps_local, sigma_n_local)
norm_res0 = jnp.linalg.norm(res)
def cond_fun(state):
norm_res, niter, _ = state
return jnp.logical_and(norm_res / norm_res0 > tol, niter < Nitermax)
def body_fun(state):
norm_res, niter, history = state
y_local, deps_local, sigma_n_local, res = history
j = drdy(y_local, deps_local, sigma_n_local)
j_inv_vp = jnp.linalg.solve(j, -res)
y_local = y_local + j_inv_vp
res = r(y_local, deps_local, sigma_n_local)
norm_res = jnp.linalg.norm(res)
history = y_local, deps_local, sigma_n_local, res
niter += 1
return (norm_res, niter, history)
history = (y_local, deps_local, sigma_n_local, res)
norm_res, niter_total, y_local = jax.lax.while_loop(cond_fun, body_fun, (norm_res0, niter, history))
sigma_local = y_local[0][:stress_dim]
dlambda = y_local[0][-1]
sigma_elas_local = C_elas @ deps_local
yielding = f(sigma_n_local + sigma_elas_local)
return sigma_local, (sigma_local, niter_total, yielding, norm_res, dlambda)
Consistent tangent stiffness matrix#
Not only is the automatic differentiation able to compute the derivative of a mathematical expression but also a numerical algorithm. For instance, AD can calculate the derivative of the function performing return-mapping with respect to its output, the stress tensor \(\boldsymbol{\sigma}\). In the context of the consistent tangent moduli \(\boldsymbol{C}_\text{tang}\), this feature becomes very useful, as there is no need to write an additional program computing the stress derivative.
JAX’s AD tool permits taking the derivative of the function return_mapping,
which is factually the while loop. The derivative is taken with respect to the
first output and the remaining outputs are used as auxiliary data. Thus, the
derivative dsigma_ddeps returns both values of the consistent tangent moduli
and the stress tensor, so there is no need in a supplementary computation of the
stress tensor.
dsigma_ddeps = jax.jacfwd(return_mapping, has_aux=True)
Defining external operator#
Once we define the function dsigma_ddeps, which evaluates both the
external operator and its derivative locally, we can simply vectorize it and
define the final implementation of the external operator derivative.
Note
The function dsigma_ddeps containing a while_loop is designed to be called
at a single Gauss point that’s why we need to vectorize it for the all points
of our functional space S. For this purpose we use the vmap function of JAX.
It creates another while_loop, which terminates only when all mapped loops
terminate. Find further details in this
discussion.
dsigma_ddeps_vec = jax.jit(jax.vmap(dsigma_ddeps, in_axes=(0, 0)))
def C_tang_impl(deps):
deps_ = deps.reshape((-1, stress_dim))
sigma_n_ = sigma_n.x.array.reshape((-1, stress_dim))
(C_tang_global, state) = dsigma_ddeps_vec(deps_, sigma_n_)
sigma_global, niter, yielding, norm_res, _dlambda = state
unique_iters, counts = jnp.unique(niter, return_counts=True)
if MPI.COMM_WORLD.rank == 0:
print("\tInner Newton summary:")
print(f"\t\tUnique number of iterations: {unique_iters}")
print(f"\t\tCounts of unique number of iterations: {counts}")
print(f"\t\tMaximum f: {jnp.max(yielding)}")
print(f"\t\tMaximum residual: {jnp.max(norm_res)}")
return C_tang_global.reshape(-1), sigma_global.reshape(-1)
Similarly to the von Mises example, we do not implement explicitly the evaluation of the external operator. Instead, we obtain its values during the evaluation of its derivative and then update the values of the operator in the main Newton loop.
def sigma_external(derivatives):
if derivatives == (1,):
return C_tang_impl
else:
raise NotImplementedError(f"No external function is defined for the requested derivative {derivatives}.")
sigma.external_function = sigma_external
Defining the forms#
q = fem.Constant(domain, default_scalar_type((0, -gamma)))
def F_ext(v):
return ufl.dot(q, v) * dx
u_hat = ufl.TrialFunction(V)
F = ufl.inner(epsilon(v), sigma) * dx - F_ext(v)
J = ufl.derivative(F, Du, u_hat)
F_replaced, F_external_operators = replace_external_operators(F)
J_replaced, J_external_operators = replace_external_operators(J)
F_form = fem.form(F_replaced)
J_form = fem.form(J_replaced)
Variables initialization and compilation#
Before solving the problem we have to initialize values of the stiffness matrix, as it requires for the system assembling. During the first loading step, we expect an elastic response only, so it’s enough to solve the constitutive equations for a relatively small displacement field at each Gauss point. This results in initializing the consistent tangent moduli with elastic ones.
Du.x.array[:] = 1.0
sigma_n.x.array[:] = 0.0
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -2.2109628558533108
Maximum residual: 0.0
Solving the problem#
Similarly to the von Mises tutorial, we use NonlinearProblem to solve the
global problem with SNES. To enable the external operators update at each
iteration of SNES before the vector and matrix assembly, we wrote a simple
wrapper assemble_residual_with_callback.
petsc_options = {
"snes_type": "vinewtonrsls",
"snes_linesearch_type": "basic",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_atol": 1.0e-8,
"snes_rtol": 1.0e-8,
"snes_max_it": 100,
"snes_monitor": "",
}
problem = NonlinearProblem(
F_replaced, Du, J=J_replaced, bcs=bcs, petsc_options_prefix="demo_mohr-coulomb_", petsc_options=petsc_options
)
def constitutive_update(
F_external_operators: list[FEMExternalOperator],
J_external_operators: list[FEMExternalOperator],
):
"""Update the constitutive model by evaluating the external operators."""
evaluated_operands = evaluate_operands(F_external_operators)
# Call `C_tang_impl`
((_, sigma_new),) = evaluate_external_operators(J_external_operators, evaluated_operands)
# Direct access to the external operator values
sigma.ref_coefficient.x.array[:] = sigma_new
assemble_residual_with_callback_ = partial(
assemble_residual_with_callback,
problem.u,
problem.F,
problem.J,
bcs,
constitutive_update, # external callback with respect to SNES
[F_external_operators, J_external_operators], # input arguments of the callback
)
problem.solver.setFunction(assemble_residual_with_callback_, problem.b)
After definition of the nonlinear problem and the Newton solver, we are ready to get the final result.
load_steps_1 = np.linspace(2, 22.9, 50)
load_steps_2 = np.array([22.96, 22.99])
load_steps = np.concatenate([load_steps_1, load_steps_2])
num_increments = len(load_steps)
results = np.zeros((num_increments + 1, 2))
x_point = np.array([[0, H, 0]])
cells, points_on_process = find_cell_by_point(domain, x_point)
for i, load in enumerate(load_steps):
q.value = load * np.array([0, -gamma])
if MPI.COMM_WORLD.rank == 0:
print(f"Load increment #{i}, load: {load}")
problem.solve()
u.x.petsc_vec.axpy(1.0, Du.x.petsc_vec)
u.x.scatter_forward()
sigma_n.x.array[:] = sigma.ref_coefficient.x.array
if len(points_on_process) > 0:
results[i + 1, :] = (-u.eval(points_on_process, cells)[0], load)
print(f"Slope stability factor: {-q.value[-1] * H / c}")
Load increment #0, load: 2.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -2.2109628558533108
Maximum residual: 0.0
0 SNES Function norm 1.195147420875e+05
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7544796613404463
Maximum residual: 0.0
1 SNES Function norm 4.925639610277e-10
Load increment #1, load: 2.426530612244898
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.1914324475785798
Maximum residual: 5.037132294934698e-10
0 SNES Function norm 4.756853159952e-02
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.8092635981099829
Maximum residual: 2.309724214262959e-16
1 SNES Function norm 3.033810662273e-02
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.35079912068642916
Maximum residual: 2.5286532809710436e-16
2 SNES Function norm 2.517197711850e-15
Load increment #2, load: 2.853061224489796
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3749 1]
Maximum f: 0.061159281086040185
Maximum residual: 1.369693540865897e-15
0 SNES Function norm 9.132432678327e-04
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3749 1]
Maximum f: 0.0923086218462843
Maximum residual: 8.86365051711617e-14
1 SNES Function norm 8.054297485286e-06
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3749 1]
Maximum f: 0.09244639358766138
Maximum residual: 9.008504901364507e-14
2 SNES Function norm 3.456153422835e-10
Load increment #3, load: 3.279591836734694
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3749 1]
Maximum f: 0.44844216521744684
Maximum residual: 1.10713013798582e-13
0 SNES Function norm 5.392112070247e-03
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.6457397002297012
Maximum residual: 5.6451657460124585e-09
1 SNES Function norm 9.222044508532e-04
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.6671698219110782
Maximum residual: 5.567802086334493e-09
2 SNES Function norm 2.306874769960e-06
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.6671780558585318
Maximum residual: 5.564613259183216e-09
3 SNES Function norm 8.267702328832e-12
Load increment #4, load: 3.706122448979592
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.6729159475629154
Maximum residual: 1.961913967726829e-09
0 SNES Function norm 7.708642137450e-03
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.8182731983589577
Maximum residual: 3.051809572597835e-09
1 SNES Function norm 1.306670391138e-04
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.8184517040414176
Maximum residual: 3.2657518352744414e-09
2 SNES Function norm 1.677202588038e-08
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 0.8184520904141586
Maximum residual: 3.2657981145976312e-09
3 SNES Function norm 2.733355749946e-15
Load increment #5, load: 4.13265306122449
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3747 3]
Maximum f: 0.8189820439289814
Maximum residual: 1.7557431539233226e-09
0 SNES Function norm 6.653604188183e-03
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3747 3]
Maximum f: 0.9125840840028121
Maximum residual: 2.890915312554209e-09
1 SNES Function norm 2.564375080374e-05
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3747 3]
Maximum f: 0.9124710788798258
Maximum residual: 2.9132791867771335e-09
2 SNES Function norm 1.145551699266e-09
Load increment #6, load: 4.559183673469388
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3744 5 1]
Maximum f: 0.9120013616400886
Maximum residual: 2.0946148405455534e-09
0 SNES Function norm 4.373771847428e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3744 5 1]
Maximum f: 1.0834716869126582
Maximum residual: 7.097372761180092e-09
1 SNES Function norm 9.084544888105e-05
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3744 5 1]
Maximum f: 1.0853634288842122
Maximum residual: 7.093969421360825e-09
2 SNES Function norm 1.607752520050e-08
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3744 5 1]
Maximum f: 1.0853639290780843
Maximum residual: 7.0939540969818396e-09
3 SNES Function norm 2.938849016761e-15
Load increment #7, load: 4.985714285714286
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3743 6 1]
Maximum f: 1.0817281941295946
Maximum residual: 6.913824904476724e-09
0 SNES Function norm 7.149688794666e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3743 5 2]
Maximum f: 1.3124494265785445
Maximum residual: 3.805485383617689e-09
1 SNES Function norm 2.058508967714e-04
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3743 5 2]
Maximum f: 1.3140149963446643
Maximum residual: 3.945091795514232e-09
2 SNES Function norm 7.527865371016e-08
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3743 5 2]
Maximum f: 1.3140158007718132
Maximum residual: 3.94511825818433e-09
3 SNES Function norm 1.318527466544e-14
Load increment #8, load: 5.412244897959184
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3740 8 1 1]
Maximum f: 1.3088422717205206
Maximum residual: 7.666114759355831e-09
0 SNES Function norm 7.696175706972e-03
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3739 8 2 1]
Maximum f: 1.5216973741715996
Maximum residual: 1.4413461891580881e-08
1 SNES Function norm 5.990449202790e-04
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3739 7 3 1]
Maximum f: 1.529233472068261
Maximum residual: 1.4268738974035388e-08
2 SNES Function norm 4.936076056872e-07
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3739 7 3 1]
Maximum f: 1.5292375934479279
Maximum residual: 1.426866345743804e-08
3 SNES Function norm 4.633461702536e-13
Load increment #9, load: 5.838775510204082
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3736 13 1]
Maximum f: 1.5243000402676858
Maximum residual: 7.688162207605254e-09
0 SNES Function norm 9.589791815981e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3735 11 4]
Maximum f: 1.7689821236208974
Maximum residual: 9.505123548537365e-09
1 SNES Function norm 3.862729558146e-04
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3735 11 4]
Maximum f: 1.7742284095842646
Maximum residual: 1.0807696643360129e-08
2 SNES Function norm 4.301313419793e-07
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3735 11 4]
Maximum f: 1.7742337488040278
Maximum residual: 1.0809538804017366e-08
3 SNES Function norm 5.348154010404e-13
Load increment #10, load: 6.265306122448979
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3731 16 3]
Maximum f: 1.7693631620667785
Maximum residual: 1.0598586281137948e-08
0 SNES Function norm 1.143203715428e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3730 15 5]
Maximum f: 1.9314433466529128
Maximum residual: 1.482621447998773e-08
1 SNES Function norm 1.584942427511e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3730 14 6]
Maximum f: 1.945981281529415
Maximum residual: 1.4853122404122691e-08
2 SNES Function norm 3.025107791995e-06
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3730 14 6]
Maximum f: 1.9460026016663963
Maximum residual: 1.486022444398606e-08
3 SNES Function norm 1.580317847763e-11
Load increment #11, load: 6.691836734693877
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3724 22 2 2]
Maximum f: 1.942723296295474
Maximum residual: 1.107082789684012e-08
0 SNES Function norm 1.107822857802e-02
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3723 18 6 3]
Maximum f: 2.122392164283419
Maximum residual: 1.7660578591258534e-08
1 SNES Function norm 5.672791184133e-04
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3723 18 6 3]
Maximum f: 2.127991178540578
Maximum residual: 1.804684437696687e-08
2 SNES Function norm 8.341174882304e-07
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3723 18 6 3]
Maximum f: 2.1279993282057252
Maximum residual: 1.804726418109548e-08
3 SNES Function norm 1.851594954092e-12
Load increment #12, load: 7.118367346938775
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3719 27 4]
Maximum f: 2.125608307652835
Maximum residual: 1.2834182578998858e-08
0 SNES Function norm 1.211225124343e-02
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3716 24 9 1]
Maximum f: 2.3289082734777344
Maximum residual: 2.829692239056086e-08
1 SNES Function norm 1.578810096150e-03
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3715 24 9 2]
Maximum f: 2.3415246200100603
Maximum residual: 1.8299013803377476e-08
2 SNES Function norm 1.892066413846e-04
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3715 24 9 2]
Maximum f: 2.3424482873879966
Maximum residual: 1.831880696662043e-08
3 SNES Function norm 3.305365333299e-08
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3715 24 9 2]
Maximum f: 2.3424484018720757
Maximum residual: 1.8318809216548028e-08
4 SNES Function norm 3.485503404078e-15
Load increment #13, load: 7.544897959183673
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3711 34 5]
Maximum f: 2.3401627518884465
Maximum residual: 2.591648055750279e-08
0 SNES Function norm 1.196820342248e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3710 28 12]
Maximum f: 2.5259332779377925
Maximum residual: 2.2868489975631276e-08
1 SNES Function norm 2.743993581195e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3710 27 13]
Maximum f: 2.5422960269574273
Maximum residual: 2.1308370545350976e-08
2 SNES Function norm 7.305829738915e-06
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3710 27 13]
Maximum f: 2.542333776222678
Maximum residual: 2.13102295844718e-08
3 SNES Function norm 7.631493803205e-11
Load increment #14, load: 7.971428571428571
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3705 42 2 1]
Maximum f: 2.5402012017918945
Maximum residual: 1.4293196563044155e-08
0 SNES Function norm 1.079588899770e-02
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3704 34 11 1]
Maximum f: 2.728556224371412
Maximum residual: 2.2875959367727578e-08
1 SNES Function norm 1.040283872579e-03
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3704 34 11 1]
Maximum f: 2.7363082848685374
Maximum residual: 2.384212431010753e-08
2 SNES Function norm 1.400314054691e-06
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3704 34 11 1]
Maximum f: 2.7363174188307755
Maximum residual: 2.3843232196162056e-08
3 SNES Function norm 2.772531732804e-12
Load increment #15, load: 8.397959183673468
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3699 49 2]
Maximum f: 2.7346869522892274
Maximum residual: 7.159754746416231e-09
0 SNES Function norm 1.091345771982e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3698 41 11]
Maximum f: 2.897365494402917
Maximum residual: 1.9947798301173393e-08
1 SNES Function norm 7.408397972280e-04
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3698 41 11]
Maximum f: 2.902394333239569
Maximum residual: 2.063408192760055e-08
2 SNES Function norm 7.993577407414e-07
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3698 41 11]
Maximum f: 2.902398634799822
Maximum residual: 2.063461756569111e-08
3 SNES Function norm 8.010292605261e-13
Load increment #16, load: 8.824489795918367
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3693 1 54 2]
Maximum f: 2.9013013029667234
Maximum residual: 7.08231966742579e-09
0 SNES Function norm 1.131724631643e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3692 53 5]
Maximum f: 3.0562073612070715
Maximum residual: 1.959059276293663e-08
1 SNES Function norm 1.008710940736e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3692 51 7]
Maximum f: 3.061188684009417
Maximum residual: 2.0280904623705412e-08
2 SNES Function norm 1.574785763099e-06
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3692 51 7]
Maximum f: 3.0611926039834336
Maximum residual: 2.0281446916682324e-08
3 SNES Function norm 3.672555190175e-12
Load increment #17, load: 9.251020408163264
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3688 61 1]
Maximum f: 3.0606359799274006
Maximum residual: 4.5392864770790375e-09
0 SNES Function norm 1.039626377063e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3687 57 6]
Maximum f: 3.193257638158007
Maximum residual: 1.383441705744067e-08
1 SNES Function norm 5.727495097146e-04
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3687 57 6]
Maximum f: 3.1964061857745647
Maximum residual: 1.4192616992203891e-08
2 SNES Function norm 4.430536443903e-07
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3687 57 6]
Maximum f: 3.196407709812544
Maximum residual: 1.4192788115922737e-08
3 SNES Function norm 2.390548387373e-13
Load increment #18, load: 9.677551020408163
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3682 67 1]
Maximum f: 3.1962792345430944
Maximum residual: 3.128096379782618e-09
0 SNES Function norm 1.077705186543e-02
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3682 63 5]
Maximum f: 3.3273846790635546
Maximum residual: 3.006978568945216e-08
1 SNES Function norm 2.242198150757e-04
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3682 63 5]
Maximum f: 3.3291427526630764
Maximum residual: 3.1822782555302277e-08
2 SNES Function norm 8.066505840979e-08
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3682 63 5]
Maximum f: 3.3291431795177773
Maximum residual: 3.182323302493181e-08
3 SNES Function norm 1.184102491236e-14
Load increment #19, load: 10.10408163265306
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3677 3 68 1 1]
Maximum f: 3.3293810538315527
Maximum residual: 2.4015415153513826e-08
0 SNES Function norm 1.009301544127e-02
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3676 69 4 1]
Maximum f: 3.446488877319546
Maximum residual: 2.190375146987787e-08
1 SNES Function norm 2.451203768063e-04
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3676 69 4 1]
Maximum f: 3.4472552731639996
Maximum residual: 1.8753968631850924e-08
2 SNES Function norm 6.772470633953e-08
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3676 69 4 1]
Maximum f: 3.447255817784177
Maximum residual: 1.8754301848669528e-08
3 SNES Function norm 1.172917492513e-14
Load increment #20, load: 10.530612244897958
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3673 75 2]
Maximum f: 3.4477457741927275
Maximum residual: 1.2300926117364232e-08
0 SNES Function norm 8.769338038872e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3672 74 4]
Maximum f: 3.553612234013521
Maximum residual: 2.0731766840377245e-08
1 SNES Function norm 1.398871589577e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3672 74 4]
Maximum f: 3.5575500941879583
Maximum residual: 2.0880012096417512e-08
2 SNES Function norm 1.379070835927e-06
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3672 74 4]
Maximum f: 3.5575533740577785
Maximum residual: 2.088004709092809e-08
3 SNES Function norm 2.265880711664e-12
Load increment #21, load: 10.957142857142856
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3669 79 2]
Maximum f: 3.5582224238556823
Maximum residual: 1.3134271505929055e-09
0 SNES Function norm 8.691186079801e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3669 2 75 4]
Maximum f: 3.6575939425921224
Maximum residual: 3.5954580547797606e-08
1 SNES Function norm 3.411985921926e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3669 1 76 4]
Maximum f: 3.660080990277343
Maximum residual: 2.9021506361322593e-09
2 SNES Function norm 2.666288752862e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3669 1 76 4]
Maximum f: 3.6600810303141427
Maximum residual: 2.902156029806634e-09
3 SNES Function norm 1.419463352103e-11
Load increment #22, load: 11.383673469387753
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3664 3 82 1]
Maximum f: 3.660889726729081
Maximum residual: 2.63979344706999e-08
0 SNES Function norm 8.784785795509e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3662 2 84 2]
Maximum f: 3.7516486387692427
Maximum residual: 4.6779433017254816e-08
1 SNES Function norm 2.334921659320e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3662 1 85 2]
Maximum f: 3.7550666405164477
Maximum residual: 4.7243832184345176e-08
2 SNES Function norm 9.229833616549e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3662 1 85 2]
Maximum f: 3.7550698944933534
Maximum residual: 4.7242393308053703e-08
3 SNES Function norm 1.329652289464e-10
Load increment #23, load: 11.810204081632651
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3660 4 84 2]
Maximum f: 3.7559609444110227
Maximum residual: 3.3906255711259566e-08
0 SNES Function norm 7.703871014509e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3659 3 85 3]
Maximum f: 3.8382982030471484
Maximum residual: 2.1412612538567907e-08
1 SNES Function norm 7.732693973153e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3659 2 86 3]
Maximum f: 3.839826825119879
Maximum residual: 2.1816109592044653e-08
2 SNES Function norm 4.461456405308e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3659 2 86 3]
Maximum f: 3.8398290986266823
Maximum residual: 2.1815665630053643e-08
3 SNES Function norm 6.876437231580e-11
Load increment #24, load: 12.23673469387755
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3653 3 91 3]
Maximum f: 3.8407568214012886
Maximum residual: 4.205786872143268e-08
0 SNES Function norm 1.146900093023e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3652 1 94 3]
Maximum f: 3.9276871850440984
Maximum residual: 5.796850157100406e-08
1 SNES Function norm 2.495998296202e-03
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3653 94 3]
Maximum f: 3.931199398279331
Maximum residual: 5.724728191399184e-08
2 SNES Function norm 1.623104017199e-05
Inner Newton summary:
Unique number of iterations: [1 3 4]
Counts of unique number of iterations: [3653 94 3]
Maximum f: 3.9312052106144804
Maximum residual: 5.724503023849471e-08
3 SNES Function norm 3.768453910687e-10
Load increment #25, load: 12.663265306122447
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3644 4 99 2 1]
Maximum f: 3.9321439660188786
Maximum residual: 5.062351073379173e-09
0 SNES Function norm 8.353184330108e-03
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3646 101 2 1]
Maximum f: 4.327729961142127
Maximum residual: 3.2333895366297856e-08
1 SNES Function norm 8.334150687846e-03
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3646 101 2 1]
Maximum f: 4.311288264135284
Maximum residual: 3.239739585803084e-08
2 SNES Function norm 4.637932352969e-06
Inner Newton summary:
Unique number of iterations: [1 3 4 5]
Counts of unique number of iterations: [3646 101 2 1]
Maximum f: 4.311303131281363
Maximum residual: 3.239702983691758e-08
3 SNES Function norm 2.609546429694e-11
Load increment #26, load: 13.089795918367345
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3639 1 108 2]
Maximum f: 4.313176238016535
Maximum residual: 2.9576101537911422e-08
0 SNES Function norm 1.194001516472e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3637 2 108 3]
Maximum f: 4.867253928816307
Maximum residual: 1.7686694220956025e-08
1 SNES Function norm 6.098877544132e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3638 2 107 3]
Maximum f: 4.959680353880248
Maximum residual: 1.6621645116617333e-08
2 SNES Function norm 1.108269377350e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3638 2 107 3]
Maximum f: 4.956079107449424
Maximum residual: 1.660921547763378e-08
3 SNES Function norm 3.073549374671e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3638 2 107 3]
Maximum f: 4.956077715187972
Maximum residual: 1.6609216731140572e-08
4 SNES Function norm 1.135287372370e-13
Load increment #27, load: 13.516326530612243
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3628 5 116 1]
Maximum f: 4.957344021789497
Maximum residual: 1.8117081073778124e-08
0 SNES Function norm 1.240366718545e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3628 3 116 3]
Maximum f: 5.683701149471084
Maximum residual: 1.29582788583561e-08
1 SNES Function norm 2.161926288287e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3627 4 115 4]
Maximum f: 5.720335933291622
Maximum residual: 1.2955192092424684e-08
2 SNES Function norm 1.004507493488e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3627 4 115 4]
Maximum f: 5.72249493226828
Maximum residual: 1.2857682905584956e-08
3 SNES Function norm 3.137931206994e-08
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3627 4 115 4]
Maximum f: 5.722495071385579
Maximum residual: 1.2857665104816402e-08
4 SNES Function norm 4.748744333229e-15
Load increment #28, load: 13.942857142857141
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3620 3 127]
Maximum f: 5.723571310910744
Maximum residual: 8.619153454244503e-09
0 SNES Function norm 1.195591532823e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3616 2 129 3]
Maximum f: 6.645080068812957
Maximum residual: 1.738732535103612e-08
1 SNES Function norm 4.874621421367e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3616 1 130 3]
Maximum f: 6.811772121665673
Maximum residual: 1.8920523131957322e-08
2 SNES Function norm 6.244285417515e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3616 1 130 3]
Maximum f: 6.812691928883815
Maximum residual: 1.892775917619562e-08
3 SNES Function norm 8.457653563827e-09
Load increment #29, load: 14.36938775510204
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3608 4 138]
Maximum f: 6.813577582394888
Maximum residual: 2.639051194295415e-08
0 SNES Function norm 1.319402690197e-02
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3607 2 141]
Maximum f: 7.596225865539337
Maximum residual: 2.7316580197235044e-08
1 SNES Function norm 6.000358681716e-03
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3607 2 141]
Maximum f: 7.662614238371001
Maximum residual: 3.175527417152072e-08
2 SNES Function norm 2.231686172695e-05
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3607 2 141]
Maximum f: 7.662669562454697
Maximum residual: 3.1765978091119795e-08
3 SNES Function norm 8.253914813013e-10
Load increment #30, load: 14.795918367346937
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3595 6 149]
Maximum f: 7.663139576851215
Maximum residual: 2.236738324419593e-08
0 SNES Function norm 1.503558807032e-02
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3594 3 153]
Maximum f: 8.454917434393016
Maximum residual: 2.0862186964093345e-08
1 SNES Function norm 6.182460517858e-03
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3594 3 153]
Maximum f: 8.561832222793038
Maximum residual: 2.6117716805807332e-08
2 SNES Function norm 5.709638472817e-05
Inner Newton summary:
Unique number of iterations: [1 2 3]
Counts of unique number of iterations: [3594 3 153]
Maximum f: 8.562225052308554
Maximum residual: 2.6136949223244306e-08
3 SNES Function norm 6.605695739847e-09
Load increment #31, load: 15.222448979591835
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3582 6 161 1]
Maximum f: 8.562438155256029
Maximum residual: 5.366138021577231e-08
0 SNES Function norm 1.356272783579e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3578 7 162 2 1]
Maximum f: 9.462706125067813
Maximum residual: 2.9615101542016138e-08
1 SNES Function norm 1.108959098908e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3578 6 163 2 1]
Maximum f: 9.531525798162972
Maximum residual: 3.6424349006276095e-08
2 SNES Function norm 2.776717567051e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3578 6 163 2 1]
Maximum f: 9.531403525835671
Maximum residual: 3.667264602785827e-08
3 SNES Function norm 1.547137372137e-09
Load increment #32, load: 15.648979591836733
Inner Newton summary:
Unique number of iterations: [1 2 3 5]
Counts of unique number of iterations: [3558 6 185 1]
Maximum f: 9.531507663600546
Maximum residual: 2.2943694420017607e-08
0 SNES Function norm 1.690023171389e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3555 2 188 4 1]
Maximum f: 11.208475801670893
Maximum residual: 1.3925568425773714e-08
1 SNES Function norm 1.136486356817e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3555 1 189 4 1]
Maximum f: 11.328808304184575
Maximum residual: 2.0711834741775247e-08
2 SNES Function norm 4.043685669749e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3555 1 189 4 1]
Maximum f: 11.328996289219967
Maximum residual: 2.0738730317501987e-08
3 SNES Function norm 3.107153995133e-09
Load increment #33, load: 16.07551020408163
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3542 2 204 2]
Maximum f: 11.329238008815347
Maximum residual: 1.1473379324306277e-08
0 SNES Function norm 1.445355802854e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3540 4 198 8]
Maximum f: 12.390016864905684
Maximum residual: 3.481616574678234e-08
1 SNES Function norm 8.639360130456e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3540 6 196 8]
Maximum f: 12.473904687662001
Maximum residual: 5.455887360562095e-08
2 SNES Function norm 1.258589840385e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3540 6 196 8]
Maximum f: 12.473963542525727
Maximum residual: 5.4478640655470284e-08
3 SNES Function norm 2.119210658480e-10
Load increment #34, load: 16.502040816326527
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3522 7 219 2]
Maximum f: 12.474040223165485
Maximum residual: 3.085274150679181e-08
0 SNES Function norm 1.409596865654e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3522 10 210 8]
Maximum f: 13.621030460090028
Maximum residual: 2.6838150764623378e-08
1 SNES Function norm 3.751124853547e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3520 10 212 8]
Maximum f: 13.672897613110022
Maximum residual: 3.5242905330548976e-08
2 SNES Function norm 1.294558689652e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3521 9 212 8]
Maximum f: 13.683042197974698
Maximum residual: 3.0360429816504824e-08
3 SNES Function norm 5.415828097352e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3521 9 212 8]
Maximum f: 13.683257893898398
Maximum residual: 3.0314245577901104e-08
4 SNES Function norm 1.090266487866e-09
Load increment #35, load: 16.928571428571427
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3507 13 229 1]
Maximum f: 13.683249738187433
Maximum residual: 1.8274113931619184e-08
0 SNES Function norm 1.773988907407e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3499 4 238 9]
Maximum f: 15.181050137105986
Maximum residual: 5.335973348691197e-08
1 SNES Function norm 4.041295762343e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3499 4 238 9]
Maximum f: 15.293687390994899
Maximum residual: 2.3828154536703193e-08
2 SNES Function norm 2.084752657070e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3499 4 238 9]
Maximum f: 15.293641206172365
Maximum residual: 2.383595316363325e-08
3 SNES Function norm 5.396666457432e-10
Load increment #36, load: 17.355102040816327
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3480 11 258 1]
Maximum f: 15.293621570476796
Maximum residual: 4.3733262304470264e-08
0 SNES Function norm 1.733879266120e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3478 7 254 11]
Maximum f: 16.658296344583988
Maximum residual: 4.011336971532695e-08
1 SNES Function norm 1.499749094517e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3478 5 256 11]
Maximum f: 16.753380479727067
Maximum residual: 1.1986831605564918e-08
2 SNES Function norm 3.007459439834e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3478 5 256 11]
Maximum f: 16.753739177855
Maximum residual: 1.1976581529300277e-08
3 SNES Function norm 1.425970722451e-09
Load increment #37, load: 17.781632653061223
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3464 11 273 1 1]
Maximum f: 16.753695524790956
Maximum residual: 5.451325624478899e-08
0 SNES Function norm 1.861829821566e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3461 11 271 5 2]
Maximum f: 18.04763005348131
Maximum residual: 6.100411597474193e-08
1 SNES Function norm 2.630504081379e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3460 8 275 5 2]
Maximum f: 18.109439698744698
Maximum residual: 7.997603669665567e-08
2 SNES Function norm 3.124568213974e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3460 8 275 5 2]
Maximum f: 18.10995427013044
Maximum residual: 8.032064151817757e-08
3 SNES Function norm 2.300727934798e-09
Load increment #38, load: 18.20816326530612
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3435 8 305 2]
Maximum f: 18.109866894786556
Maximum residual: 3.5438443056582133e-08
0 SNES Function norm 1.800669384053e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3427 8 309 6]
Maximum f: 19.682636417527156
Maximum residual: 5.462722615856093e-08
1 SNES Function norm 4.274620709122e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3427 7 310 6]
Maximum f: 19.816059584913678
Maximum residual: 7.360482487709092e-08
2 SNES Function norm 2.487642596250e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3427 7 310 6]
Maximum f: 19.816699769298737
Maximum residual: 7.313228879539957e-08
3 SNES Function norm 1.096984795880e-09
Load increment #39, load: 18.63469387755102
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3407 15 325 3]
Maximum f: 19.816608603382683
Maximum residual: 3.294772775452879e-08
0 SNES Function norm 1.931304808040e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3406 7 328 9]
Maximum f: 21.389004861044647
Maximum residual: 4.9589345093117305e-08
1 SNES Function norm 5.092741146638e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3406 7 328 9]
Maximum f: 21.48012370240415
Maximum residual: 5.53855898231622e-08
2 SNES Function norm 2.203884275220e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3406 7 328 9]
Maximum f: 21.480128176158512
Maximum residual: 5.5395368107855424e-08
3 SNES Function norm 1.122075644037e-09
Load increment #40, load: 19.061224489795915
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3385 13 350 2]
Maximum f: 21.48001805230716
Maximum residual: 3.935054878743307e-08
0 SNES Function norm 1.807185057934e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3377 14 351 8]
Maximum f: 23.101864004083005
Maximum residual: 7.489324996645121e-08
1 SNES Function norm 4.096957992682e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3376 8 358 8]
Maximum f: 23.18123342797546
Maximum residual: 5.4376246706206966e-08
2 SNES Function norm 3.801885713101e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3376 7 359 8]
Maximum f: 23.185942545326498
Maximum residual: 5.475290252338967e-08
3 SNES Function norm 2.086410212481e-08
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3376 7 359 8]
Maximum f: 23.18594259385968
Maximum residual: 5.4752907245139965e-08
4 SNES Function norm 6.317456239223e-15
Load increment #41, load: 19.487755102040815
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3356 11 381 2]
Maximum f: 23.18582201266842
Maximum residual: 6.081588052723994e-08
0 SNES Function norm 2.318768587420e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3352 6 380 11 1]
Maximum f: 25.49716808894195
Maximum residual: 4.983964059444414e-08
1 SNES Function norm 5.854691572877e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3351 6 382 10 1]
Maximum f: 25.567903421736972
Maximum residual: 7.710589688417007e-08
2 SNES Function norm 1.857969440167e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3351 6 382 10 1]
Maximum f: 25.57006462133442
Maximum residual: 7.630891785529839e-08
3 SNES Function norm 1.324273255742e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3351 6 382 10 1]
Maximum f: 25.570065383153835
Maximum residual: 7.630882498462073e-08
4 SNES Function norm 4.503982142656e-14
Load increment #42, load: 19.91428571428571
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3324 11 413 1 1]
Maximum f: 25.56996105387789
Maximum residual: 5.606058391898205e-08
0 SNES Function norm 1.991950861898e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3320 6 413 10 1]
Maximum f: 27.485399219692063
Maximum residual: 4.0171595610476786e-08
1 SNES Function norm 6.102183392930e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3320 8 411 10 1]
Maximum f: 27.605560170484807
Maximum residual: 3.6863139067903996e-08
2 SNES Function norm 5.558752205564e-05
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3320 8 411 10 1]
Maximum f: 27.605715757299244
Maximum residual: 3.684919797754657e-08
3 SNES Function norm 7.661646727445e-09
Load increment #43, load: 20.34081632653061
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3291 14 444 1]
Maximum f: 27.605586022820038
Maximum residual: 5.77554611634033e-08
0 SNES Function norm 2.099035542902e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3288 11 443 8]
Maximum f: 29.843867326126606
Maximum residual: 7.139782457659531e-08
1 SNES Function norm 4.733736371878e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3289 11 442 8]
Maximum f: 29.969805340763685
Maximum residual: 6.912245793708935e-08
2 SNES Function norm 5.313476427743e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3289 11 442 8]
Maximum f: 29.9740324334093
Maximum residual: 6.694398006386842e-08
3 SNES Function norm 1.094543345653e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3289 11 442 8]
Maximum f: 29.97403261643221
Maximum residual: 6.694298781530843e-08
4 SNES Function norm 1.306026357878e-14
Load increment #44, load: 20.767346938775507
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3259 20 469 2]
Maximum f: 29.973884886574155
Maximum residual: 6.014568428884466e-08
0 SNES Function norm 2.246453731859e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3253 15 472 10]
Maximum f: 32.90205667394817
Maximum residual: 4.9787387694762134e-08
1 SNES Function norm 4.783821924576e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3254 9 477 10]
Maximum f: 33.07703435009244
Maximum residual: 5.3298261180874925e-08
2 SNES Function norm 5.291518412053e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3254 9 477 10]
Maximum f: 33.08102796048997
Maximum residual: 5.337579548551162e-08
3 SNES Function norm 4.903824806473e-08
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3254 9 477 10]
Maximum f: 33.0810279767684
Maximum residual: 5.337579907135085e-08
4 SNES Function norm 8.114326890943e-15
Load increment #45, load: 21.193877551020407
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3222 14 512 1 1]
Maximum f: 33.08088303623554
Maximum residual: 4.586197473990603e-08
0 SNES Function norm 2.150367183099e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3214 14 509 12 1]
Maximum f: 36.42903952633916
Maximum residual: 6.060362972975903e-08
1 SNES Function norm 8.415327791298e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3212 14 512 11 1]
Maximum f: 36.59087082197401
Maximum residual: 7.033600288801433e-08
2 SNES Function norm 2.624314374599e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3212 14 512 11 1]
Maximum f: 36.592839261624775
Maximum residual: 7.041867461508215e-08
3 SNES Function norm 2.176387425440e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3212 14 512 11 1]
Maximum f: 36.59284038499766
Maximum residual: 7.041868404508342e-08
4 SNES Function norm 1.347894271744e-13
Load increment #46, load: 21.620408163265303
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3178 21 549 1 1]
Maximum f: 36.59268395743359
Maximum residual: 2.3037654988075924e-08
0 SNES Function norm 2.628890154187e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3175 15 548 12]
Maximum f: 40.880163997826415
Maximum residual: 4.861051666770778e-08
1 SNES Function norm 9.704135791354e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3174 13 551 12]
Maximum f: 41.09688044587984
Maximum residual: 5.888597584182205e-08
2 SNES Function norm 3.940138059941e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3174 13 551 12]
Maximum f: 41.09933447019429
Maximum residual: 5.8958580774101356e-08
3 SNES Function norm 3.404476110573e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3174 13 551 12]
Maximum f: 41.09933731996514
Maximum residual: 5.8958638573954166e-08
4 SNES Function norm 3.270177049172e-13
Load increment #47, load: 22.046938775510203
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3138 21 588 3]
Maximum f: 41.09916240554589
Maximum residual: 1.5081275654740336e-07
0 SNES Function norm 2.587384322553e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3129 27 582 12]
Maximum f: 46.39966838688594
Maximum residual: 6.262153546401841e-08
1 SNES Function norm 6.838178297472e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3127 21 589 12 1]
Maximum f: 46.93592390989389
Maximum residual: 1.029980579453794e-07
2 SNES Function norm 1.362075338948e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3127 19 591 12 1]
Maximum f: 46.972122873918885
Maximum residual: 8.753638526277173e-08
3 SNES Function norm 8.572437202811e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3127 19 591 12 1]
Maximum f: 46.972210151266566
Maximum residual: 8.753996181047612e-08
4 SNES Function norm 1.261714820928e-10
Load increment #48, load: 22.4734693877551
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3086 30 627 7]
Maximum f: 46.97203139227406
Maximum residual: 6.7433309935667e-08
0 SNES Function norm 2.944549916639e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3072 26 633 19]
Maximum f: 54.79142280645419
Maximum residual: 1.2807527127899936e-07
1 SNES Function norm 1.035617501699e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3072 25 634 19]
Maximum f: 55.92681732819188
Maximum residual: 1.9443039193793728e-07
2 SNES Function norm 1.664047906161e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3072 23 637 18]
Maximum f: 55.97127030749246
Maximum residual: 1.9966760045480525e-07
3 SNES Function norm 6.382437507972e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [3072 23 637 18]
Maximum f: 55.97133470774807
Maximum residual: 1.996764497174674e-07
4 SNES Function norm 1.046913089773e-10
Load increment #49, load: 22.9
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [3006 36 697 10 1]
Maximum f: 55.97111640595372
Maximum residual: 9.181422348345212e-08
0 SNES Function norm 3.451445654693e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2979 20 730 20 1]
Maximum f: 71.24902596407796
Maximum residual: 1.7736360700474487e-07
1 SNES Function norm 3.791796477487e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2971 20 726 32 1]
Maximum f: 78.51083340525379
Maximum residual: 3.422498489782812e-07
2 SNES Function norm 6.547900362150e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2973 18 725 33 1]
Maximum f: 79.67788407053496
Maximum residual: 5.530154988050549e-07
3 SNES Function norm 1.411257682785e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2973 18 725 33 1]
Maximum f: 79.7736403579503
Maximum residual: 5.717286832215061e-07
4 SNES Function norm 1.114569889288e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2973 18 725 33 1]
Maximum f: 79.77518101584543
Maximum residual: 5.719374306199863e-07
5 SNES Function norm 1.307185596190e-08
Inner Newton summary:
Unique number of iterations: [1 2 3 4 5]
Counts of unique number of iterations: [2973 18 725 33 1]
Maximum f: 79.77518100703828
Maximum residual: 5.719374229450116e-07
6 SNES Function norm 1.488788036129e-14
Load increment #50, load: 22.96
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2875 27 823 25]
Maximum f: 79.77483948791325
Maximum residual: 7.72099727319028e-08
0 SNES Function norm 4.510532966358e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2980 163 590 17]
Maximum f: 18.040683536843076
Maximum residual: 3.2623013119252314e-08
1 SNES Function norm 3.813905199850e-01
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2956 143 646 5]
Maximum f: 21.192127984270947
Maximum residual: 3.221171524363646e-08
2 SNES Function norm 2.731881915507e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2939 126 676 9]
Maximum f: 27.358774625347547
Maximum residual: 3.970759463402072e-08
3 SNES Function norm 9.518819283938e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2937 103 700 10]
Maximum f: 32.93695984912313
Maximum residual: 7.979458984561214e-08
4 SNES Function norm 2.440368481888e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2935 101 703 11]
Maximum f: 34.048736088512754
Maximum residual: 2.4912423837556754e-08
5 SNES Function norm 4.704398167596e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2935 101 703 11]
Maximum f: 34.14222362720718
Maximum residual: 2.4406672810829435e-08
6 SNES Function norm 1.117573929995e-06
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2935 101 703 11]
Maximum f: 34.14235387239141
Maximum residual: 2.440610684836912e-08
7 SNES Function norm 5.486436202416e-12
Load increment #51, load: 22.99
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2923 112 712 3]
Maximum f: 34.14219556832545
Maximum residual: 6.430010385698741e-08
0 SNES Function norm 1.458357152583e-02
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2925 161 662 2]
Maximum f: 29.88026032182711
Maximum residual: 2.9275667379816684e-08
1 SNES Function norm 5.932855542944e-03
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2926 157 665 2]
Maximum f: 31.32107837332672
Maximum residual: 3.334337741320375e-08
2 SNES Function norm 6.233209758949e-04
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2926 157 665 2]
Maximum f: 31.363803938673307
Maximum residual: 3.3202357084262685e-08
3 SNES Function norm 4.410931139292e-07
Inner Newton summary:
Unique number of iterations: [1 2 3 4]
Counts of unique number of iterations: [2926 157 665 2]
Maximum f: 31.363876920814583
Maximum residual: 3.3202640421228784e-08
4 SNES Function norm 1.805132554050e-12
Slope stability factor: 6.663768115942029
Verification#
Critical load#
According to Chen and Liu [1990], we can derive analytically the slope stability factor \(l_\text{lim}\) for the standard Mohr-Coulomb plasticity model (without apex smoothing) under plane strain assumption for associative plastic flow
\[
l_\text{lim} = \gamma_\text{lim} H/c,
\]
where \(\gamma_\text{lim}\) is an associated value of the soil self-weight. In particular, for the rectangular slope with the friction angle \(\phi\) equal to \(30^\circ\), \(l_\text{lim} = 6.69\) [Chen and Liu, 1990]. Thus, by computing \(\gamma_\text{lim}\) from the formula above, we can progressively increase the second component of the gravitational body force \(\boldsymbol{q}=[0, -\gamma]^T\), up to the critical value \(\gamma_\text{lim}^\text{num}\), when the perfect plasticity plateau is reached on the loading-displacement curve at the \((0, H)\) point and then compare \(\gamma_\text{lim}^\text{num}\) against analytical \(\gamma_\text{lim}\).
By demonstrating the loading-displacement curve on the figure below we approve that the yield strength limit reached for \(\gamma_\text{lim}^\text{num}\) is close to \(\gamma_\text{lim}\).
if len(points_on_process) > 0:
l_lim = 6.69
gamma_lim = l_lim / H * c
plt.plot(results[:, 0], results[:, 1], "o-", label=r"$\gamma$")
plt.axhline(y=gamma_lim, color="r", linestyle="--", label=r"$\gamma_\text{lim}$")
plt.xlabel(r"Displacement of the slope $u_x$ at $(0, H)$ [mm]")
plt.ylabel(r"Soil self-weight $\gamma$ [MPa/mm$^3$]")
plt.grid()
plt.legend()
The slope profile reaching its stability limit:
try:
import pyvista
print(pyvista.global_theme.jupyter_backend)
import dolfinx.plot
W = fem.functionspace(domain, ("Lagrange", 1, (gdim,)))
u_tmp = fem.Function(W, name="Displacement")
u_tmp.interpolate(u)
plotter = pyvista.Plotter(window_size=[600, 400], off_screen=True)
topology, cell_types, x = dolfinx.plot.vtk_mesh(domain)
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
vals = np.zeros((x.shape[0], 3))
vals[:, : len(u_tmp)] = u_tmp.x.array.reshape((x.shape[0], len(u_tmp)))
grid["u"] = vals
warped = grid.warp_by_vector("u", factor=20)
plotter.add_mesh(warped, show_edges=False, show_scalar_bar=False)
plotter.view_xy()
plotter.camera.tight()
image = plotter.screenshot(None, transparent_background=True, return_img=True)
plt.imshow(image)
plt.axis("off")
except ImportError:
print("pyvista required for this plot")
static
2026-07-14 19:08:33.513 ( 0.288s) [ 7FE996E21140]vtkXOpenGLRenderWindow.:1460 WARN| bad X server connection. DISPLAY=:99.0
Yield surface#
We verify that the constitutive model is correctly implemented by tracing the yield surface. We generate several stress paths and check whether they remain within the Mohr-Coulomb yield surface. The stress tracing is performed in the Haigh-Westergaard coordinates \((\xi, \rho, \theta)\) which are defined as follows
\[
\xi = \frac{1}{\sqrt{3}}I_1, \quad \rho =
\sqrt{2J_2}, \quad \sin(3\theta) = -\frac{3\sqrt{3}}{2}
\frac{J_3}{J_2^{3/2}},
\]
where \(J_3(\boldsymbol{\sigma}) = \det(\boldsymbol{s})\) is the third invariant of the deviatoric part of the stress tensor, \(\xi\) is the deviatoric coordinate, \(\rho\) is the radial coordinate and the angle \(\theta \in [-\frac{\pi}{6}, \frac{\pi}{6}]\) is called Lode or stress angle.
To generate the stress paths we use the principal stresses formula written in Haigh-Westergaard coordinates as follows
\[\begin{split}
\begin{pmatrix}
\sigma_{I} \\
\sigma_{II} \\
\sigma_{III}
\end{pmatrix}
= p
\begin{pmatrix}
1 \\
1 \\
1
\end{pmatrix}
+ \frac{\rho}{\sqrt{2}}
\begin{pmatrix}
\cos\theta + \frac{\sin\theta}{\sqrt{3}} \\
-\frac{2\sin\theta}{\sqrt{3}} \\
\frac{\sin\theta}{\sqrt{3}} - \cos\theta
\end{pmatrix},
\end{split}\]
where \(p = \xi/\sqrt{3}\) is a hydrostatic variable and \(\sigma_{I} \geq \sigma_{II} \geq \sigma_{III}\).
Now we generate the loading path by evaluating principal stresses in Haigh-Westergaard coordinates for the Lode angle \(\theta\) being varied from \(-\frac{\pi}{6}\) to \(\frac{\pi}{6}\) with fixed \(\rho\) and \(p\).
N_angles = 50
N_loads = 9 # number of loadings or paths
eps = 0.00001
R = 0.7 # fix the values of rho
p = 0.1 # fix the deviatoric coordinate
theta_1 = -np.pi / 6
theta_2 = np.pi / 6
theta_values = np.linspace(theta_1 + eps, theta_2 - eps, N_angles)
theta_returned = np.empty((N_loads, N_angles))
rho_returned = np.empty((N_loads, N_angles))
sigma_returned = np.empty((N_loads, N_angles, stress_dim))
# fix an increment of the stress path
dsigma_path = np.zeros((N_angles, stress_dim))
dsigma_path[:, 0] = (R / np.sqrt(2)) * (np.cos(theta_values) + np.sin(theta_values) / np.sqrt(3))
dsigma_path[:, 1] = (R / np.sqrt(2)) * (-2 * np.sin(theta_values) / np.sqrt(3))
dsigma_path[:, 2] = (R / np.sqrt(2)) * (np.sin(theta_values) / np.sqrt(3) - np.cos(theta_values))
sigma_n_local = np.zeros_like(dsigma_path)
sigma_n_local[:, 0] = p
sigma_n_local[:, 1] = p
sigma_n_local[:, 2] = p
derviatoric_axis = tr
Then, we define and vectorize functions rho, Lode_angle and sigma_tracing
evaluating respectively the coordinates \(\rho\), \(\theta\) and the corrected (or
“returned”) stress tensor for a certain stress state. sigma_tracing calls the
function return_mapping, where the constitutive model was defined via JAX
previously.
def rho(sigma_local):
s = dev @ sigma_local
return jnp.sqrt(2.0 * J2(s))
def Lode_angle(sigma_local):
s = dev @ sigma_local
arg = -(3.0 * jnp.sqrt(3.0) * J3(s)) / (2.0 * jnp.sqrt(J2(s) * J2(s) * J2(s)))
arg = jnp.clip(arg, -1.0, 1.0)
angle = 1.0 / 3.0 * jnp.arcsin(arg)
return angle
def sigma_tracing(sigma_local, sigma_n_local):
deps_elas = S_elas @ sigma_local
sigma_corrected, state = return_mapping(deps_elas, sigma_n_local)
yielding = state[2]
return sigma_corrected, yielding
Lode_angle_v = jax.jit(jax.vmap(Lode_angle, in_axes=(0)))
rho_v = jax.jit(jax.vmap(rho, in_axes=(0)))
sigma_tracing_v = jax.jit(jax.vmap(sigma_tracing, in_axes=(0, 0)))
For each stress path, we call the function sigma_tracing_v to get the
corrected stress state and then we project it onto the deviatoric plane \((\rho,
\theta)\) with a fixed value of \(p\).
for i in range(N_loads):
print(f"Loading path#{i}")
dsigma, yielding = sigma_tracing_v(dsigma_path, sigma_n_local)
dp = dsigma @ tr / 3.0 - p
dsigma -= np.outer(dp, derviatoric_axis) # projection on the same deviatoric plane
sigma_returned[i, :] = dsigma
theta_returned[i, :] = Lode_angle_v(dsigma)
rho_returned[i, :] = rho_v(dsigma)
print(f"max f: {jnp.max(yielding)}\n")
sigma_n_local[:] = dsigma
Loading path#0
max f: -2.005661796528811
Loading path#1
max f: -1.647414005283729
Loading path#2
max f: -1.208026577509537
Loading path#3
max f: -0.7355419142072797
Loading path#4
max f: -0.24734105482689195
Loading path#5
max f: 0.24936204365846093
Loading path#6
max f: 0.5729074243253178
Loading path#7
max f: 0.6673099640326825
Loading path#8
max f: 0.6947128480833951
Then, by knowing the expression of the standrad Mohr-Coulomb yield surface in principle stresses, we can obtain an analogue expression in Haigh-Westergaard coordinates, which leads us to the following equation:
(7)#\[
\frac{\rho}{\sqrt{6}}(\sqrt{3}\cos\theta + \sin\phi
\sin\theta) - p\sin\phi - c\cos\phi= 0.
\]
Thus, we restore the standard Mohr-Coulomb yield surface:
def MC_yield_surface(theta_, p):
"""Restores the coordinate `rho` satisfying the standard Mohr-Coulomb yield
criterion."""
rho = (np.sqrt(2) * (c * np.cos(phi) + p * np.sin(phi))) / (
np.cos(theta_) - np.sin(phi) * np.sin(theta_) / np.sqrt(3)
)
return rho
rho_standard_MC = MC_yield_surface(theta_values, p)
Finally, we plot the yield surface:
colormap = cm.plasma
colors = colormap(np.linspace(0.0, 1.0, N_loads))
fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(8, 8))
# Mohr-Coulomb yield surface with apex smoothing
for i, color in enumerate(colors):
rho_total = np.array([])
theta_total = np.array([])
for j in range(12):
angles = j * np.pi / 3 - j % 2 * theta_returned[i] + (1 - j % 2) * theta_returned[i]
theta_total = np.concatenate([theta_total, angles])
rho_total = np.concatenate([rho_total, rho_returned[i]])
ax.plot(theta_total, rho_total, ".", color=color)
# standard Mohr-Coulomb yield surface
theta_standard_MC_total = np.array([])
rho_standard_MC_total = np.array([])
for j in range(12):
angles = j * np.pi / 3 - j % 2 * theta_values + (1 - j % 2) * theta_values
theta_standard_MC_total = np.concatenate([theta_standard_MC_total, angles])
rho_standard_MC_total = np.concatenate([rho_standard_MC_total, rho_standard_MC])
ax.plot(theta_standard_MC_total, rho_standard_MC_total, "-", color="black")
ax.set_yticklabels([])
norm = mcolors.Normalize(vmin=0.1, vmax=0.7 * 9)
sm = plt.cm.ScalarMappable(cmap=colormap, norm=norm)
sm.set_array([])
cbar = fig.colorbar(sm, ax=ax, orientation="vertical")
cbar.set_label(r"Magnitude of the stress path deviator, $\rho$ [MPa]")
plt.show()
Each colour represents one loading path. The circles are associated with the loading during the elastic phase. Once the loading reaches the elastic limit, the circles start outlining the yield surface, which in the limit lay along the standard Mohr-Coulomb one without smoothing (black contour).
Taylor test#
Here, we perform a Taylor test to check that the form \(F\) and its Jacobian \(J\)
are consistent zeroth- and first-order approximations of the residual \(F\). In
particular, the test verifies that the program dsigma_ddeps_vec
obtained by the JAX’s AD returns correct values of the external operator
\(\boldsymbol{\sigma}\) and its derivative \(\boldsymbol{C}_\text{tang}\), which
define \(F\) and \(J\) respectively.
To perform the test, we introduce the operators \(\mathcal{F}: V \rightarrow V^\prime\) and \(\mathcal{J}: V \rightarrow \mathcal{L}(V, V^\prime)\) defined as follows:
\[
\langle \mathcal{F}(\boldsymbol{u}), \boldsymbol{v} \rangle :=
F(\boldsymbol{u}; \boldsymbol{v}), \quad \forall \boldsymbol{v} \in V,
\]
\[
\langle (\mathcal{J}(\boldsymbol{u}))(k\boldsymbol{\delta u}),
\boldsymbol{v} \rangle := J(\boldsymbol{u}; k\boldsymbol{\delta u},
\boldsymbol{v}), \quad \forall \boldsymbol{v} \in V,
\]
where \(V^\prime\) is a dual space of \(V\), \(\langle \cdot, \cdot \rangle\) is the \(V^\prime \times V\) duality pairing and \(\mathcal{L}(V, V^\prime)\) is a space of bounded linear operators from \(V\) to its dual.
Then, by following the Taylor’s theorem on Banach spaces and perturbating the functional \(\mathcal{F}\) in the direction \(k \, \boldsymbol{δu} \in V\) for \(k > 0\), the zeroth and first order Taylor reminders \(r_k^0\) and \(r_k^1\) have the following mesh-independent convergence rates in the dual space \(V^\prime\):
(8)#\[
\| r_k^0 \|_{V^\prime} := \| \mathcal{F}(\boldsymbol{u} + k \,
\boldsymbol{\delta u}) - \mathcal{F}(\boldsymbol{u}) \|_{V^\prime}
\longrightarrow 0 \text{ at } O(k),
\]
(9)#\[
\| r_k^1 \|_{V^\prime} := \| \mathcal{F}(\boldsymbol{u} + k \,
\boldsymbol{\delta u}) - \mathcal{F}(\boldsymbol{u}) - \,
(\mathcal{J}(\boldsymbol{u}))(k\boldsymbol{\delta u}) \|_{V^\prime}
\longrightarrow 0 \text{ at } O(k^2).
\]
In order to compute the norm of an element \(f \in V^\prime\) from the dual space \(V^\prime\), we apply the Riesz representation theorem, which states that there is a linear isometric isomorphism \(\mathcal{R} : V^\prime \to V\), which associates a linear functional \(f\) with a unique element \(\mathcal{R} f = \boldsymbol{u} \in V\). In practice, within a finite subspace \(V_h \subset V\), the Riesz map \(\mathcal{R}\) is represented by the matrix \(\mathsf{L}^{-1}\), the inverse of the Laplacian operator [Kirby, 2010]
\[
\mathsf{L}_{ij} = \int\limits_\Omega \nabla\varphi_i \cdot \nabla\varphi_j \mathrm{d} x , \quad i,j = 1, \dots, n,
\]
where \(\{\varphi_i\}_{i=1}^{\dim V_h}\) is a set of basis function of the space \(V_h\).
If the Euclidean vectors \(\mathsf{r}_k^i \in \mathbb{R}^{\dim V_h}, \, i \in \{0,1\}\) represent the Taylor remainders from (8)–(9) in the finite space, then the dual norms are computed through the following formula [Kirby, 2010]
(10)#\[
\| r_k^i \|^2_{V^\prime_h} = (\mathsf{r}_k^i)^T \mathsf{L}^{-1} \mathsf{r}_k^i, \quad i \in \{0,1\}.
\]
In practice, the vectors \(\mathsf{r}_k^i\) are defined through the residual vector \(\mathsf{F} \in \mathbb{R}^{\dim V_h}\) and the Jacobian matrix \(\mathsf{J} \in \mathbb{R}^{\dim V_h\times\dim V_h}\)
(11)#\[
\mathsf{r}_k^0 = \mathsf{F}(\mathsf{u} + k \, \mathsf{\delta u}) - \mathsf{F}(\mathsf{u}) \in \mathbb{R}^n,
\]
(12)#\[
\mathsf{r}_k^1 = \mathsf{F}(\mathsf{u} + k \, \mathsf{\delta u}) -
\mathsf{F}(\mathsf{u}) - \, \mathsf{J}(\mathsf{u}) \cdot k\mathsf{\delta
u} \in \mathbb{R}^n,
\]
where \(\mathsf{u} \in \mathbb{R}^{\dim V_h}\) and \(\mathsf{\delta u} \in \mathbb{R}^{\dim V_h}\) represent dispacement fields \(\boldsymbol{u} \in V_h\) and \(\boldsymbol{\delta u} \in V_h\).
Now we can proceed with the Taylor test implementation. Let us first start with defining the Laplace operator.
L_form = fem.form(ufl.inner(ufl.grad(u_hat), ufl.grad(v)) * ufl.dx)
L = fem.petsc.assemble_matrix(L_form, bcs=bcs)
L.assemble()
Riesz_solver = PETSc.KSP().create(domain.comm)
Riesz_solver.setType("preonly")
Riesz_solver.getPC().setType("lu")
Riesz_solver.setOperators(L)
y = fem.Function(V, name="Riesz_representer_of_r") # r - a Taylor remainder
Now we initialize main variables of the plasticity problem.
# Reset main variables to zero including the external operators values
sigma_n.x.array[:] = 0.0
sigma.ref_coefficient.x.array[:] = 0.0
J_external_operators[0].ref_coefficient.x.array[:] = 0.0
# Reset the values of the consistent tangent matrix to elastic moduli
Du.x.array[:] = 1.0
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -2.2109628558533108
Maximum residual: 0.0
As the derivatives of the constitutive model are different for elastic and plastic phases, we must consider two initial states for the Taylor test. For this reason, we solve the problem once for a certain loading value to get the initial state close to the one with plastic deformations but still remain in the elastic phase.
i = 0
load = 2.0
q.value = load * np.array([0, -gamma])
Du.x.array[:] = 1e-8
if MPI.COMM_WORLD.rank == 0:
print(f"Load increment #{i}, load: {load}")
problem.solve()
u.x.petsc_vec.axpy(1.0, Du.x.petsc_vec)
u.x.scatter_forward()
sigma_n.x.array[:] = sigma.ref_coefficient.x.array
Du0 = np.copy(Du.x.array)
sigma_n0 = np.copy(sigma_n.x.array)
Load increment #0, load: 2.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -2.210962855861672
Maximum residual: 0.0
0 SNES Function norm 5.506494297258e-02
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7544796600182155
Maximum residual: 0.0
1 SNES Function norm 2.938490495611e-14
If we take into account the initial stress state sigma_n0 computed in the cell
above, we perform the Taylor test for the plastic phase, otherwise we stay in
the elastic one.
Finally, we define the function perform_Taylor_test, which returns the norms
of the Taylor reminders in dual space (10)–(12).
k_list = np.logspace(-2.0, -6.0, 5)[::-1]
def perform_Taylor_test(Du0, sigma_n0):
# r0 = F(Du0 + k*δu) - F(Du0)
# r1 = F(Du0 + k*δu) - F(Du0) - k*J(Du0)*δu
Du.x.array[:] = Du0
sigma_n.x.array[:] = sigma_n0
evaluated_operands = evaluate_operands(F_external_operators)
((_, sigma_new),) = evaluate_external_operators(J_external_operators, evaluated_operands)
sigma.ref_coefficient.x.array[:] = sigma_new
F0 = fem.petsc.assemble_vector(F_form) # F(Du0)
F0.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(F0, bcs)
J0 = fem.petsc.assemble_matrix(J_form, bcs=bcs)
J0.assemble() # J(Du0)
Ju = J0.createVecLeft() # Ju = J0 @ u
δu = fem.Function(V)
δu.x.array[:] = Du0 # δu == Du0
zero_order_remainder = np.zeros_like(k_list)
first_order_remainder = np.zeros_like(k_list)
for i, k in enumerate(k_list):
Du.x.array[:] = Du0 + k * δu.x.array
evaluated_operands = evaluate_operands(F_external_operators)
((_, sigma_new),) = evaluate_external_operators(J_external_operators, evaluated_operands)
sigma.ref_coefficient.x.array[:] = sigma_new
F_delta = fem.petsc.assemble_vector(F_form) # F(Du0 + h*δu)
F_delta.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(F_delta, bcs)
J0.mult(δu.x.petsc_vec, Ju) # Ju = J(Du0)*δu
Ju.scale(k) # Ju = k*Ju
r0 = F_delta - F0
r1 = F_delta - F0 - Ju
Riesz_solver.solve(r0, y.x.petsc_vec) # y = L^{-1} r0
y.x.scatter_forward()
zero_order_remainder[i] = np.sqrt(r0.dot(y.x.petsc_vec)) # sqrt{r0^T L^{-1} r0}
Riesz_solver.solve(r1, y.x.petsc_vec) # y = L^{-1} r1
y.x.scatter_forward()
first_order_remainder[i] = np.sqrt(r1.dot(y.x.petsc_vec)) # sqrt{r1^T L^{-1} r1}
return zero_order_remainder, first_order_remainder
print("Elastic phase")
zero_order_remainder_elastic, first_order_remainder_elastic = perform_Taylor_test(Du0, 0.0)
print("Plastic phase")
zero_order_remainder_plastic, first_order_remainder_plastic = perform_Taylor_test(Du0, sigma_n0)
Elastic phase
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7544796600182155
Maximum residual: 0.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7544777931845448
Maximum residual: 0.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7544609916686915
Maximum residual: 0.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7542929752409684
Maximum residual: 0.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7526126841444585
Maximum residual: 0.0
Inner Newton summary:
Unique number of iterations: [1]
Counts of unique number of iterations: [3750]
Maximum f: -0.7357971890466226
Maximum residual: 0.0
Plastic phase
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.1914324503846303
Maximum residual: 5.037132209749364e-10
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.191434439616415
Maximum residual: 5.037159079560956e-10
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.1914523427045762
Maximum residual: 5.037359395191501e-10
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.1916313737948552
Maximum residual: 5.039397845836618e-10
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.1934217055527925
Maximum residual: 5.05980367704671e-10
Inner Newton summary:
Unique number of iterations: [1 3]
Counts of unique number of iterations: [3748 2]
Maximum f: 1.211327098969814
Maximum residual: 5.266651628781508e-10
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].loglog(k_list, zero_order_remainder_elastic, "o-", label=r"$\|r_k^0\|_{V^\prime}$")
axs[0].loglog(k_list, first_order_remainder_elastic, "o-", label=r"$\|r_k^1\|_{V^\prime}$")
annotation.slope_marker((2e-4, 5e-5), 1, ax=axs[0], poly_kwargs={"facecolor": "tab:blue"})
axs[0].text(0.5, -0.2, "(a) Elastic phase", transform=axs[0].transAxes, ha="center", va="top")
axs[1].loglog(k_list, zero_order_remainder_plastic, "o-", label=r"$\|r_k^0\|_{V^\prime}$")
annotation.slope_marker((2e-4, 5e-5), 1, ax=axs[1], poly_kwargs={"facecolor": "tab:blue"})
axs[1].loglog(k_list, first_order_remainder_plastic, "o-", label=r"$\|r_k^1\|_{V^\prime}$")
annotation.slope_marker((2e-4, 5e-13), 2, ax=axs[1], poly_kwargs={"facecolor": "tab:orange"})
axs[1].text(0.5, -0.2, "(b) Plastic phase", transform=axs[1].transAxes, ha="center", va="top")
for i in range(2):
axs[i].set_xlabel("k")
axs[i].set_ylabel("Taylor remainder norm")
axs[i].legend()
axs[i].grid()
plt.tight_layout()
plt.show()
first_order_rate = np.polyfit(np.log(k_list), np.log(zero_order_remainder_elastic), 1)[0]
second_order_rate = np.polyfit(np.log(k_list), np.log(first_order_remainder_elastic), 1)[0]
print(f"Elastic phase:\n\tthe 1st order rate = {first_order_rate:.2f}\n\tthe 2nd order rate = {second_order_rate:.2f}")
first_order_rate = np.polyfit(np.log(k_list), np.log(zero_order_remainder_plastic), 1)[0]
second_order_rate = np.polyfit(np.log(k_list[1:]), np.log(first_order_remainder_plastic[1:]), 1)[0]
print(f"Plastic phase:\n\tthe 1st order rate = {first_order_rate:.2f}\n\tthe 2nd order rate = {second_order_rate:.2f}")
Elastic phase:
the 1st order rate = 1.00
the 2nd order rate = 0.00
Plastic phase:
the 1st order rate = 1.00
the 2nd order rate = 1.93
For the elastic phase (a) the zeroth-order Taylor remainder \(r_k^0\) achieves the first-order convergence rate, whereas the first-order remainder \(r_k^1\) is computed at the level of machine precision due to the constant Jacobian. Similarly to the elastic flow, the zeroth-order Taylor remainder \(r_k^0\) of the plastic phase (b) reaches the first-order convergence, whereas the first-order remainder \(r_k^1\) achieves the second-order convergence rate, as expected.