Add option for resetting failing Newton solver

[finsberg]

Sometimes when you solve a nonlinear problem, the solver might diverge which leaves the solver in a bad state.

Consider the following scenario

solver = NewtonSolver(...)
target_value = 1.0
step = 0.1
while parameter.value < target_value:
    parameter.value += step
    try:
        solver.solve()
    except RuntimeError:
        # Do smaller step and reset solver
        parameter.value -= step
        step /= 2.0

However, when the solver fails it might be in a bad state which needs to be reset to the last good state. My suggestion would be to allow for resetting the state in a seamless way, e.g

solver = NewtonSolver(...)
target_value = 1.0
step = 0.1
while parameter.value < target_value:
    good_state = solver.state
    parameter.value += step
    try:
        solver.solve()
    except RuntimeError:
        parameter.value -= step
        step /= 2.0
        solver.reset_state(good_state)

or even implement a method called try_to_solve_or_reset (or just make it default to reset the state).

[jorgensd]

Can’t you simply control this with adding some checks, and change error_on_no convergence to False?

a set_state (rather than a reset state) could be ok with the appropriate documentation (that data is copied out of an array, rather than re-assigned in the ufl forms.

[finsberg]

My current workaround is as follows

states = [u, p]
solver = scifem.NewtonSolver(F, J, states, bcs=bcs)

old_states = []
for f in states:
    f_old = dolfinx.fem.Function(u.function_space)
    f_old.x.array[:] = f.x.array
    old_states.append(f_old)

target_value = 1.0
step = 0.1
while parameter.value < target_value:
    parameter.value += step
    try:
        solver.solve()
    except RuntimeError:
        parameter.value -= step
        step /= 2.0
        for f, f_old in zip(states, old_states):
            f.x.array[:] = f_old.x.array
      
        # Create a new solver
        solver = scifem.NewtonSolver(F, J, states, bcs=bcs)
    else:
        for f, f_old in zip(states, old_states):
            f_old.x.array[:] = f.x.array

I think it would be nice if we could do this without the need for creating a new solver.

[jorgensd]

Maybe I am missing the point, but what exactly has to be reset (except the unknown ?)