Scaling may not reach Residual in time

[ibaned]

In reasoning about #76, @bgranzow and I are concerned that if the Residual is requested before the Jacobian, it will not receive the same scaling as the Jacobian. This issue tracks my progress investigating that. If it is true, it would give incorrect results with Scaling on. @ikalash

[ibaned]

Some initial debugging suggests that our concern is valid:

DEBUG: computing residual in step 0, scaling it by scaleVec from step -2
DEBUG: computing residual in step 0, scaling it by scaleVec from step -2
DEBUG: computing residual in step 1, scaling it by scaleVec from step -2
DEBUG: computing jacobian in step 1, scaling it by scaleVec from step 1
DEBUG: computing residual in step 1, scaling it by scaleVec from step 1
DEBUG: computing residual in step 1, scaling it by scaleVec from step 1
DEBUG: computing residual in step 2, scaling it by scaleVec from step 1
DEBUG: computing jacobian in step 2, scaling it by scaleVec from step 2
DEBUG: computing residual in step 2, scaling it by scaleVec from step 2
DEBUG: computing residual in step 2, scaling it by scaleVec from step 2
DEBUG: computing residual in step 3, scaling it by scaleVec from step 2
DEBUG: computing jacobian in step 3, scaling it by scaleVec from step 3
DEBUG: computing residual in step 3, scaling it by scaleVec from step 3
DEBUG: computing residual in step 3, scaling it by scaleVec from step 3
DEBUG: computing residual in step 4, scaling it by scaleVec from step 3
DEBUG: computing jacobian in step 4, scaling it by scaleVec from step 4
DEBUG: computing residual in step 4, scaling it by scaleVec from step 4
DEBUG: computing residual in step 4, scaling it by scaleVec from step 4
DEBUG: computing residual in step 5, scaling it by scaleVec from step 4
DEBUG: computing jacobian in step 5, scaling it by scaleVec from step 5
DEBUG: computing residual in step 5, scaling it by scaleVec from step 5
DEBUG: computing residual in step 5, scaling it by scaleVec from step 5

Here, "step" is a monotonic numbering of all nonlinear solver steps in the simulation. Step -2 represents the identity scaling vector that is used if the Jacobian has not yet been computed. The problem is very simple: the residual is requested before the Jacobian, but the Jacobian is needed to compute the scaling for the Residual. Thus, each residual gets the scaling that is appropriate to the Jacobian from the previous nonlinear step, or unity if there is no previous step. This output is using the Line Search Based nonlinear solver, which computes the residual twice after solving the linear system to check convergence. Unfortunately, the first residual evaluation is the critical one that is used to form the linear system:

Start of Continuation Step 2 : Parameter: Time = 2.000e+00 from 1.000e+00
Continuation Method: Natural
Current step size  = 1.000e+00   Previous step size = 1.000e+00
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

DEBUG: computing residual in step 2, scaling it by scaleVec from step 1

************************************************************************
-- Nonlinear Solver Step 0 --
||F|| = 8.096e+00  step = 0.000e+00  dx = 0.000e+00
************************************************************************

DEBUG: computing jacobian in step 2, scaling it by scaleVec from step 2

  *******************************************************
  ***** Belos Iterative Solver:  Block Gmres
  ***** Maximum Iterations: 200
  ***** Block Size: 1
  ***** Residual Test:
  *****   Test 1 : Belos::StatusTestImpResNorm<>: (2-Norm Res Vec) / (2-Norm Res0), tol = 1e-06
  *******************************************************
  Iter   0, [ 1] :    1.000000e+00
  Iter   1, [ 1] :    3.114679e-02
  Iter   2, [ 1] :    2.908969e-03
  Iter   3, [ 1] :    2.721727e-04
  Iter   4, [ 1] :    2.910549e-05
  Iter   5, [ 1] :    2.843900e-06
  Iter   6, [ 1] :    3.113096e-07
DEBUG: computing residual in step 2, scaling it by scaleVec from step 2

************************************************************************
-- Nonlinear Solver Step 1 --
||F|| = 5.470e-06  step = 1.000e+00  dx = 8.876e+00 (Converged!)
************************************************************************

************************************************************************
-- Final Status Test Results --
Converged....OR Combination ->
  **...........Model Evaluator Flag
  Converged....F-Norm = 5.470e-06 < 1.000e-03
               (Unscaled Two-Norm, Absolute Tolerance)
************************************************************************

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
End of Continuation Step 2 : Parameter: Time = 2.000e+00 from 1.000e+00
--> Step Converged in 1 Nonlinear Solver Iterations!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

DEBUG: computing residual in step 2, scaling it by scaleVec from step 2

  Calling Predictor with method: Constant

[rppawlo]

We are very careful about this in the nox solvers and make sure to recalculate norms (instead of reusing across Newton iterations) since weights can change even between Newton steps. Are you using the nox PrePostOperator to update the weights? See, for example, the file:

Trilinos/packages/nox/src-thyra/NOX_PrePostOperator_RowSumScaling.H

This object allows you to specify when to update the scaling so that it can be used for both steady-state (pre-iterate) and transient (pre-solve) solves.

[ibaned]

@rppawlo Right now all the scaling logic is hardcoded deep in Albany, and the scaling happens before we return the Residual and Jacobian from our Thyra ModelEvaluator. The scaling factors are computed when the Jacobian is computed, at every Newton step. I'd certainly like to use an existing Trilinos system, because Scaling should really not be the concern of the physics application. PrePostOperators look promising. Is there an example of activating this through NOX ParameterLists ?