Newton residual and tangent normalization for multiple unknowns
Interacting with the non-linear solver when we are dealing with different unknowns that have different amplitudes (for instance a displacement-pressure formulation, where displacements could be of the order of 10^{-3}
and pressures of order 10^4
) must be done with care.
Indeed, the convergence criteria for the solver are of three different kinds:
- The norm of the residual is smaller than a fixed value (absolute convergence)
- The norm of the residual is smaller than
N
times the residual of the first Newton iteration (relative convergence) - The norm of the increment of the solution is smaller than a fixed value (absolute step size convergence)
One of the criteria could be satisfied prior to reaching the real solution of the problem.
In order to circumvent this, we need to rescale the residual and the tangent by the magnitude value of each unknown and by the magnitude of the residual itself.
Typically for a displacement/pressure formulation we would solve:
\begin{pmatrix}
\frac{\partial R_u^{k}}{\partial u^{k}} & \frac{\partial R_u^{k}}{\partial p^{k}} & \\
\frac{\partial R_p^{k}}{\partial u^{k}} & \frac{\partial R_p^{k}}{\partial p^{k}} &
\end{pmatrix}
\begin{pmatrix}
\Delta u^{k+1} & \\
\Delta p^{k+1} &
\end{pmatrix}
=
\begin{pmatrix}
R_u^{k} & \\
R_p^{k} &
\end{pmatrix},
where k
is the previous Newton iteration and k+1
the current iteration.
Note that in reality we would have displacement and pressure nodes intertwinned to limit the bandwidth of the tangent matrix.
Then the solution is vector is updated as
\begin{pmatrix}
u^{k+1} & \\
p^{k+1} &
\end{pmatrix}
=
\begin{pmatrix}
\Delta u^{k+1} & \\
\Delta p^{k+1} &
\end{pmatrix}
+
\begin{pmatrix}
u^{k} & \\
p^{k} &
\end{pmatrix}.
The rescaling that we need to do is the following:
Solve the system
Ax = B
where we have
A = \begin{pmatrix}
1 / |R_u^{0}|& 0 & \\
0 & 1 / |R_p^{0}| &
\end{pmatrix}
\begin{pmatrix}
\frac{\partial R_u^{k}}{\partial u^{k}} & \frac{\partial R_u^{k}}{\partial p^{k}} & \\
\frac{\partial R_p^{k}}{\partial u^{k}} & \frac{\partial R_p^{k}}{\partial p^{k}} &
\end{pmatrix}
\begin{pmatrix}
u_{mag} & 0 & \\
0 & p_{mag} &
\end{pmatrix},
\quad
x =
\begin{pmatrix}
u_{mag} & 0 & \\
0 & p_{mag} &
\end{pmatrix}
\begin{pmatrix}
\Delta u^{k+1} & \\
\Delta p^{k+1} &
\end{pmatrix},
\quad
B = \begin{pmatrix}
R_u^{k} / |R_u^{0}| & \\
R_p^{k} / |R_p^{0}| &
\end{pmatrix}.
That way the convergence will be monitored properly. Note that for the division by the magnitude of the residuals |R_u^{k}|
, we will need to check that the numerical value is not too small in order to avoid numerical errors. Typically if we have a local value where NumericNS::IsZero()
returns true then we should not rescale the value in question (division by 1. would do the trick).
In short what we need to implement this is:
- a new magnitude field associated to each unknown
- rescale the non-linear system (rhs, matrix and state), which will be done through numbering subset interpolations.
The best case scenario would be that if the magnitude fields exist and that more than one unknown belongs to the numbering subset for which the system is solved then the rescaling would be done under the hood.