Linear solver residual tolerances should not always be normalized
Functionality to add/problem to solve
In the implementation of the linear solvers, the calculated L2-norm (?) of the residual vector is compared to user supplied absolute and relative tolerances. In the code, before comparison, the calculated residual is normalized with an estimate of the scale of the problem, which is computed from an initial guess for the solution and the RHS of the equation system. Actually, this turns the "absolute" tolerance into some kind of relative tolerance.
This is helpful (though misleading), if the expected absolute values of the residuals are unknown to the user. It is a disadvantage, if the expected absolute errors are known.
Thus it is proposed to optionally switch of the normalization of the residuals. A proposed patch is attached.
IMHO, the current specification of a tolerance and a relTol in the solver dictionary is weird, because actually both tolerances are compared to normalized residuals in the code. One by some properties of the equation system, the other by the initial residual. The treatment in the code "looks" different, but actually both variants are normalizations. But cleaning this up would require larger changes in the linear solvers, that I didn't want to tackle.
If the expected residuals and tolerances are well known, it is beneficial to use absolute tolerances for the solvers and not normalized ones, because it enhances the performance. Example: Transient tracer transport. In a simulation period with no or few tracer in the system, in the current implementation the normalization factor is near zero. Thus, the normalization of the residuals with "near zero" results in excessive iterations to reach the tolerance. The acquired accuracy is gained at high cost and might be unnecessary, if later in the transient development, at a much higher level of tracer concentrations, the interesting things happen.
Make the normalization of residuals a user choice. The proposed patch is an enhancement in lduMatrix::solver::normFactor . A user specified "residualNormalization" can then be added to the solver dictionary. If omitted, the old behaviour is retained. If set to "none", the normalization method returns "1", i.e. normalization is switched off. The currently implemented method is also available as axbScaling". For future use, other choices can easily be added.
What does success look like, and how can we measure that?
Problem: The necessary effort (i.e. iterations) for the linear solvers should depend on the dynamics of the physical processes. The current normalization hinders that.
Example: Unzip the example and run scalarTransportFoam. It is a simple example for a tracer transport, with a tracer injection starting at t=0.1 s. In the time < 0.1 s, only the initial value of tracer is in the system (very small) and is unchanged. Thus the solver should spend minimal effort. In the current implementation, the first initial residual is rescaled to 1 and unexplicably drops until t=0.1 s. The linear solver spends 5-8 unnecessary iterations in each time step. This is unwanted. After the tracer injection started, the number of iterations stays in the same range, while it should be higher then.
With the patch: Uncomment "residualNormalization" in system/fvSolution/solvers and re-run scalarTransportFoam. Result: The solver spends zero iterations before the injection starts, because the initial residual are very low. After injection (t=0.1s), the iterations increase.
Links / references
The requested functionality has been implemented and tested. Both are attached (see above).
The diff of lduMatrixSolver.C is here: