Hi everyone,

Currently I am trying to debug one of a simulation of a time dependent 3D heat equation on a fine cylindrical mesh. I am using very small time steps (0.01 seconds) for this problem due to the fast dynamics of the system, and whilst my code appears to run well for half of the time steps with the residual norm decreasing steadily, but then suddenly, I get a error message:

Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 87 iterations (PETSc reason DIVERGED_DTOL, residual norm ||r|| = 6.781851e+12).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 8
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown

I should note that this code compiles fully in less dynamic systems (those that change on the order of seconds).

What is causing this error and how can it be fixed?

Below are the parameters for the solver I have used for this system, which I suspect are causing the problem.

problem = NonlinearVariationalProblem(F, T, [], Gain)
	solver  = NonlinearVariationalSolver(problem)

	solver.parameters["newton_solver"]["relative_tolerance"] = 1E-5
	solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-5
	solver.parameters["newton_solver"]["convergence_criterion"] = "residual"
	solver.parameters["newton_solver"]["error_on_nonconvergence"] = True
	solver.parameters["newton_solver"]["linear_solver"] = "bicgstab"
	#solver.parameters["newton_solver"]["lu_solver"]["symmetric"] = True 
	solver.parameters["newton_solver"]["maximum_iterations"] = 100
	solver.parameters["newton_solver"]["preconditioner"] = "jacobi"
	#solver.parameters["newton_solver"]["relaxation_parameter"] = 1.05
	solver.parameters["newton_solver"]["report"] = True

	solver.parameters["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True 
	solver.parameters["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-5
	solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-5
	solver.parameters["newton_solver"]["krylov_solver"]["monitor_convergence"] = True
	solver.parameters["newton_solver"]["krylov_solver"]["maximum_iterations"] = 100000

Hi, I’m not sure about what is happening, but let me give you some possible options:

  1. The heat equations actually suffers from numerical instabilities from ‘not too low’ diffusion, so that if the time-step is too small, the reaction term (inertia) is dominant, so maybe you need an additional bubble dof for stability.
  2. The Jacobi preconditioner isn’t very scalable, try something like petsc_amg or hypre_euclid.
  3. Your krylov solver is not reading the params you are giving. In fact, it stops at 87 iterations, so maybe try setting them through the general dolfin parameters interface? I’m not sure on that one.

Hope it helps somehow. Best!

In addition to the Nicolas’ suggestions, maybe you can try using direct solver, just to see if the problem is not caused by iterative solver (although, based on your code, it is possible that you have already tried this before).

Alternatively, you can try using PETScSNESSolver(), which has a bit more options than built-in NewtonSolver, or you can designing your own NewtonSolver as proposed here.

Regarding the solver parameters, you can pass them alternatively in the following way (example is for PETScSNESSolver()):

nonlinear_solver = PETScSNESSolver()
PETScOptions.set("snes_rtol ", "1e-4")
PETScOptions.set("snes_atol ", "1e-4")
PETScOptions.set("ksp_initial_guess_nonzero", "true")
PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "hypre")

Beware that for some reason setting the relative and absolute tolerance for nonlinear solver in this way does not work, so you will need to do it as in your example. Hope this helps.