Hello !
I have a problem with some parameter n. I need get a solution for a certain value of n. I need to increase it slowly otherwise the solvers fails at some point.
So I’d like a way to increase n as quick as possible. My idea is having an increasing step and if solver fails, decrease the step and roll back to the old solution.
However I don’t know how to get back to the old state properly.
Example
Let’s take the hyperelastic demo Hyperelasticity — FEniCSx tutorial and modify the end of the script (see below). Here I make two small n increment (0 to 1 and 1 to 2), followed a big n increment (2 to 30) to get a failure. Then I try to revert to n=1 state and increment from 1 to 2 again, but the solver failed.
What am I doing wrong ?
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, nls
L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 2))
def left(x): return np.isclose(x[0], 0)
def right(x): return np.isclose(x[0], L)
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full(len(left_facets), 1, dtype=np.int32), np.full(len(right_facets), 2, dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
v = ufl.TestFunction(V)
u = fem.Function(V)
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
P = ufl.diff(psi, F)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
problem = fem.petsc.NonlinearProblem(F, u, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.error_on_nonconvergence = False # No error
tval0 = -1.5
T.value[2] = 1 * tval0
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
print ("Solved with n=1")
uSave = fem.Function(V)
uSave.x.array[:]=u.x.array[:]
T.value[2] = 2 * tval0
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
print ("Solved with n=2")
T.value[2] = 30 * tval0
num_its, converged = solver.solve(u)
assert( not converged)
u.x.array[:]=uSave.x.array[:]
print ("Back to n=1 solution")
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
print ("Solved with n=2 again")
Raphaël