Dear all,
in my time dependent problem, I am using a time step size control like in the following pseudo-code. For the sake of simplicity I decided to not submit the full running example that would reproduce the issue.
In the case of no convergence, the time step dt is reduced to a small number and the problem solved again. In the residual res as well as the tangent dres the solution w is initialized with the previously converged solution wrestart.
This works fine until the residual becomes nan in the last Newton iteration at some time step. In the output you can see that after restarting the solver, all subsequent Newton iterations become nan as well starting from the first iteration.
I am wondering whether the tangent is recomputed in the new solver call or if the previous tangent is reused for initialization.
Note: I am working on a very small time scale which could possibly be an issue. However, solving time steps with dt~4e-10 is not a problem in intermediate steps.
I am using fenicsx v0.9.0 compiled in a spack environment.
from mpi4py import MPI
from dolfinx import fem
import ufl
comm = MPI.COMM_WORLD
w = fem.Function(W)
wm1 = fem.Function(W)
wrestart = fem.Function(W)
dw = ufl.TestFunction(W)
dwt = ufl.TrialFunction(W)
delta_t = fem.Constant(domain, 5e-1)
res = something + ufl.dot( (w-wm1)/delta_t, dw )*ufl.dx # rate dependent
dres = ufl.derivative(res, w, dwt)
problem = fem.petsc.NonlinearProblem(res, w, bcs, dres)
solver = nls.petsc.NewtonSolver(comm, problem)
ksp = solver.krylov_solver
opts = PETSc.Options()
opts["nls_solve_ksp_type"] = "preonly"
opts["nls_solve_pc_type"] = "lu"
opts["nls_solve_pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
while True:
comm.barrier()
ufl.replace(res, {w:w})
ufl.replace(dres, {w:w})
ufl.replace(res, {wm1:wm1})
ufl.replace(dres, {wm1:wm1})
converged = False
restart = False
try:
(iters, converged) = solver.solve(w)
except RuntimeError: # no convergence
dt = 8e-10
restart = True
dt = comm.bcast(dt, root=0)
restart = comm.bcast(restart, root=0)
delta_t.value = dt # update dt
if restart:
t = t_restart + dt # time of previous step + new time increment
w.interpolate(wrestart)
w.x.scatter_forward()
else:
# update fields
w.x.scatter_forward()
wm1.interpolate(w)
wm1.x.scatter_forward()
# save solution
wrestart.interpolate(w)
wrestart.x.scatter_forward()
# update t
t_restart = deepcopy(t)
t += dt
# update boundary condition for next time step
...
==================================================
time step: 44
time = 2.2042e-05
with dt = 1.6384e-06
Newton iteration 0: r (abs) = 38.50881950091842 (tol = 1e-09), r (rel) = 0.2965375639781686 (tol = -inf)
PETSc Krylov solver starting to solve system.
Newton iteration 1: r (abs) = 0.006858700378961545 (tol = 1e-09), r (rel) = 5.196527462848138e-05 (tol = -inf)
...
Newton iteration 4: r (abs) = 3.2821394966861324e+16 (tol = 1e-09), r (rel) = 248672883914640.53 (tol = -inf)
...
Newton iteration 5: r (abs) = 6.787565973900231e+53 (tol = 1e-09), r (rel) = 5.14263213734473e+51 (tol = -inf)
...
Newton iteration 6: r (abs) = 1.0252211361998453e+132 (tol = 1e-09), r (rel) = 7.767637446442153e+129 (tol = -inf)
...
Newton iteration 7: r (abs) = -nan (tol = 1e-09), r (rel) = -nan (tol = -inf)
...
Newton iteration 11: r (abs) = -nan (tol = 1e-09), r (rel) = -nan (tol = -inf)
----------------
converged: False
restarting
----------------
decreasing dt = 8e-10 < 1.6384e-06
==================================================
time step: 44
time = 2.0404e-05
with dt = 8.0000e-10
Newton iteration 0: r (abs) = 0.01938575949236433 (tol = 1e-09), r (rel) = 0.00014687714293402734 (tol = -inf)
PETSc Krylov solver starting to solve
Newton iteration 1: r (abs) = -nan (tol = 1e-09), r (rel) = -nan (tol = -inf)