Solver fails, how to recover previous state

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

As you are solving a Non-linear problem, it is dependent on the initial values in u as it is used as an initial guess for the non-linear problem and the residual computation. Therefore, you would need to store the solution from the previous step, and insert it into u if you want to revert the solver to its previous state.

1 Like

But is it not what I am doing with u.x.array[:]=uSave.x.array[:] ?
In the example, I was thinking I’m doing n from 2 to 3 twice in the exact same manner, but it fails the second time. So I must be wrong somewhere.

Hm, you are right that you are resetting u.x.array.
To me it seems like there is a bug in PETSc (as I’ve check all the matrices/vectors prior to assembly, and they are the same as they were at the first instance).

I’m not quite sure how to proceed. The optimal thing would be to make a minimal example and create an issue on PETSc’s gitlab page. However, it can be quite tedious to do so. I’ll see if I can come up with any ideas as to how to simplify the problem.

1 Like

I will try to investigate this more into details the following days.
Probably, I will use the docker development environment, build dolfinx and petsc in the release with debug info mode, and use gdb to see what’s going on.
@dokken Have you send any bug report to the PETSc development community ? If not, I might do so.

Ive not had time to do so. Feel free to do it.

I am facing a similar issue of needed to revert to a previous solution (in the context of simple adaptive time stepping in elastodynamics). It looks like the only difference in our general approaches is that I am saving the checkpoint as an numpy array instead of store a full fem.Function instance. I was wondering if you had any updates on your approach.

@dokken would the quick and dirty fix be to re-instantiate the problem/solver each time it fails or does that seem impractical (re. memory, time, etc)? If I did this, would I be able to specify the initial guess as the previous save state?

Any help would be greatly appreciated!

As this is an issue with PETSc, someone would need to report it to their issue tracker.

You could re-create the solver, but you should avoid doing it too often, as shown in: Force clear PETSc solver from memory? - #8 by dokken
As the solver does not cache the previous state (it depends on the form, it should be possible to do

To give a short update, I’m now using the SNES PETSc solver instead of the basic NewtonSolver of fenicsx, and it works, so I have a feeling that the issue might not be on PETSc. Please find some code below, adapt the solver() method with your problem.

class SNESProblem():
    def __init__(self, F, u, bc):
        V = u.function_space
        du = ufl.TrialFunction(V)
        self.L = fem.form(F)
        self.a = fem.form(ufl.derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.vector)
        self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        fem.petsc.assemble_vector(F, self.L)
        fem.petsc.apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(F, [self.bc], x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        fem.petsc.assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()


def solver():
    problem = SNESProblem(F, u, bcs[0])
    b = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
    J = fem.petsc.create_matrix(problem.a)

    snes = PETSc.SNES().create()
    snes.setFunction(problem.F, b)
    snes.setJacobian(problem.J, J)
    opts = PETSc.Options()
    #opts['snes_monitor'] = None
    opts['snes_linesearch_type'] = 'basic'
    snes.setFromOptions()

    #snes.setType("newtontr")
    snes.setType("newtonls")

    max_it = 40
    snes.setTolerances(rtol=1.0e-8, max_it=max_it)
    snes.getKSP().setType("preonly")
    snes.getKSP().setTolerances(rtol=1.0e-8)
    snes.getKSP().getPC().setType("lu")
    return snes

Hi @RaphaelC! I have tried to implement your solution but I am still having the same issue. Would you be able to share the MWE from your original post with the updated solution to show the full implementation?