Implement underrelaxation

[paulromano]

Moving issue 21 from gitlab originally submitted by @aprilnovak:

Enrico currently uses standard Picard iteration, where a field in iteration i + 1 is computed as the result of some operator A acting on the previous iteration field: u(i + 1) = A(u(i), v(i)). When Doppler feedback is strong and/or the fluid flowrate is sufficiently low, oscillations might be observed between successive iterations in the heat source and/or heat transfer fields. Underrelaxation is a common approach for both reducing these oscillations and speeding up slow convergence.

Underrelaxation computes a field in iteration i + 1 as a linear combination of the previous iterate and the most recently-computed approximation as u(i + 1)=(1 - a) * u(i) + A(u(i), v(i)), where a is a constant parameter in the range (0, 1]. Because underrelaxation is so commonly-used, it would be nice to have the ability to use it in Enrico. Implementing underrelaxation involves several steps:

• Create example demonstrating oscillations (without relaxation)
• Create example demonstrating slow convergence (without relaxation)
• Implement underrelaxation capabilities. People most commonly apply underrelaxation to the heat source, so I will assume relaxation of that field, with possible future extensions to other options.
• Apply underrelaxation to the two examples to demonstrate both the correct implementation and advantages of the approach.

The above tasks will assume a simple constant-a relaxation. The next steps can be to add stochastic approximation methods, which vary a (and in some cases the number of particles simulated in each Picard iteration).

• Implement Robbins-Monro stochastic approximation (variable a, constant number of particles/iteration)
• Implement Dufek-Gudowski stochastic approximation (variable a and number of particles/iteration)
• Apply Robbins-Monro to two examples
• Apply Dufek-Gudowski to two examples