Enable iterative viscosity dampening

[naliboff]

This [WIP] PR adds an option to average (dampen) the viscosity (or other properties) between nonlinear iterations, which has been reported to help improve convergence behavior.

At present, the PR is only partially complete as the method for storing (compositional field) and updating (reaction terms) the viscosity between nonlinear iterations does not work with standard compositional fields.

Here is a summary of what I think the issue is, borrowing from a comment by @jdannberg in #4923.

Ignoring advection, the reaction term update is C_new = C_old + R. Here, R is the reaction term and C_old is the old solution, which is not updated with nonlinear iterations.

The equation I would like to implement is C_new = C_old + R = C_old + (-C_old + new_viscosity), where C_old is the viscosity is the previous nonlinear iteration and new_viscosity is the viscosity calculated in the most recent nonlinear iteration.

So, the first issue is that C_old is not updated during the nonlinear iterations. Can this be circumvented in some capacity by implementing a new custom field method?

The second issue is that the C_new values here are jumping between the same two values on subsequent nonlinear iterations, which is an issue @jdannberg also pointed out in the same comment.. However, I'm not sure why that is happening here as R is not equal to C_old

Hopefully there is path forward here.