Make the Newton method compatible with the GMG solver

[YiminJin]

There is some trouble when using the Newton method with GMG solver. The trouble was firstly discovered by @naliboff when testing #5580 with the kaus_2010 benchmark. It turns out that the modified Newton method performs poorly with GMG solver. I have run this benchmark on a AVX512 machine with 8 MPI processes and made comparisons between #5580 and the main branch. The results are:

Code version	Cell averaging on viscosity	Linear solver	Total wall time (s)	Total nonlinear iterations
Main branch	none	block AMG	811	592
Main branch	harmonic	block AMG	640	451
Main branch	harmonic	block GMG	678	460
PR5580	none	block AMG	290	171
PR5580	harmonic	block AMG	371	203
PR5580	harmonic	block GMG	719	460

The prm file is
d1e20-par-cav-test1.prm.txt
As seen, the cell averaging on viscosity has a positive effect in the main branch, while in PR5580 the situation is just opposite: both the GMG solver and the cell averaging have negative effects on the convergence rate. Notice that without cell averaging, PR5580 is the fastest and cuts off more than half of the wall time comparing with the main branch, which proves that the modified Newton solver is efficient. The question is why the modified Newton solver does not work with cell averaging and GMG method.

The reason why the cell averaging affects the Newton solver is easy to understand. When assembling the system matrix and RHS, the derivatives of viscosity are calculated IN function material_model.evaluate(), but function MaterialAveraging::average() is called AFTER function material_model.evaluate(), which changes the values of viscosity. Thus, the system matrix no longer represents the system tangent. To fix this problem, we may add a function to MaterialModel::Interface that calculates the viscosity derivatives after cell averaging, or just do something in function MaterialAveraging::average() to make the viscosity derivatives compatible with the averaged viscosity.

However, I haven't figured out why the GMG method affects the performance of the Newton solver so seriously, because I'm not familiar with the GMG method. Since the GMG method is in general faster than the block AMG method and will become the default solver of ASPECT, I think it is necessary to fix this problem. Please give me a hand if anyone has got any idea.

[YiminJin]

Thanks to @gassmoeller for merging #5580 to the main branch, which partially solves this issue. To make the GMG solver perform better with Newton method, 2 more improvements are required:

1. Implement the GMG preconditioner for the system Jacobian (it is in a TODO state now);
2. Modify the matrix-free framework so that the material average is not calculated each and every time when function local_apply() is called.

Both of these changes need a lot of work. For now I shall close this issue and wait for another chance to dive into this task.