chalmersplasmatheory/DREAM

Check magnitude of residual in Newton solver

Closed this issue · 2 comments

In a few weird cases, it is found that the Newton solver converges (in the sense that it stops updating the solution, dx -> 0) while the residual remains far from zero. These cases have inconsistent "solutions" and are even unphysical, and should be guarded against by adding an explicit check of the residual in the Newton solver.

(How to exactly check whether the residual is small remains to be figured out -- we must know the magnitude of the terms in the equation, not just the magnitude of the unknown quantity)

After checking the residual in one case, and thinking some, it seems that a feasible check could be to

  • add a special case DREAM::FVM::Operator::AddTerm(TransientTerm*), which behaves the same as the EquationTerm* counterpart, but also sets a flag indicating that a TransientTerm is present in the operator. This should also be propagated to UnknownQuantityEquation.
  • In ConvergenceChecker, add a new method for checking if each element of the residual is within tolerance. The method will need to take a pointer to the residual as input (preferably the "preconditioned" residual). The convergence checker should account for possible transient terms in the equation by dividing the relevant residual elements by the time step, dt.
  • Add a flag to the Python interface which allows enforcing a tolerance on the residual explicitly, i.e. if the residual does not meet the desired tolerance, quit. By default, only a soft check (not crashing the code) should be performed, which issues a warning if the residual is not converged.
  • Save a flag to the output which, for each time step, indicates whether the residual was adequately close to zero.

In addition to the above, the following minor fixes could also be included in the same pull request:

  • Make it possible to store the iteration progress to the output file (i.e. the text showing how different unknowns converge in each Newton iteration when setVerbose(True) has been called)
  • Add a function to py/cli/cli.py which allows one to easily translate a time in seconds to a time index, e.g. index2time().
  • A bug was discovered with how SaveDebugInfoAfter() is called. This method is called after the solution is calculated, but before it is actually updated in the solver, and so SaveDebugInfoAfter() actually works with data for the previous iteration.