Account for cases where the A block is not symmetric

[gassmoeller]

While looking into #5607 I found it weird that enabling the free surface would deteriorate the convergence behavior of the inner A block solver as much as it does. What was particularly weird was that the model would converge just fine with more cheap iterations. That made me look into the stabilization term that we apply for the free surface, namely here. This term is not necessarily symmetric with respect to the indices of the test functions. The crucial part is (scratch.phi_u[i]*g_hat) * (scratch.phi_u[j]*n_hat). You can imagine that scratch.phi_u[j] may be orthogonal to n_hat but not g_hat, which would result in different stabilization terms added if you flip the indices i and j. However, we choose CG as the solver for the A block in solver.cc, and CG requires the matrix to be symmetric. Using a non-symmetric matrix with a CG solver could be the reason for the convergence problems if the non-symmetric part of the matrix is large compared to the symmetric part.

So I tried replacing CG with FGMRES and the failing test checkpoint_07_enable_free_surface_resume would immediately pass with the solver settings that would fail with CG.

I think we can consider two options:

1. Try to detect if the A block is not symmetric and then switch to a different solver for the A block, this is what I implemented in this PR. This is the straightforward approach, but it requires us to check for all reasons why the A block can become not symmetric (e.g. what about anisotropic viscosity?). Also it seems like the non-symmetric part of the A block is typically small, otherwise we would have seen more failures of the A block solver in the past. Thus, always using FGMRES for free surface models could cost performance, while not increasing the stability for most models.
2. We could always try the CG solver and catch the exception if it fails and try again with FGMRES. This option is nice, because it simplifies the coding and we cannot accidentally 'miss' cases where the matrix becomes non-symmetric. However, it means that in some models we will try CG again for each outer GMRES iteration, even if it fails every time. Maybe we can store this state in the preconditioner (e.g. if CG fails once, always use GMRES as long as this preconditioner exists).

Has anyone else seen models where the inner A block solver would not converge in free surface models? I would be grateful for feedback and opinions.

[bangerth]

Anisotropic viscosity still leads to a symmetric matrix. You have a term of the form (eps_i, X eps_j) where X is a symmetric rank-4 tensor that must satisfy (eps_i, X eps_j) = (X eps_i, eps_j).

[bangerth]

As for the term added to the matrix, are we sure that it is actually correct? It might be worth going back to the paper in which it was proposed and double-check.

[gassmoeller]

I discussed this with @tjhei last week and I think this is the right path to go. As far as we know this is the right stabilization term unless we want to derive a completely new stabilization (which would likely require a new free-surface paper). I also implemented the switch for GMG now, and added an input parameter to force the new solver. This is useful, because we have occasionally observed models that crash in the A block preconditioner in the past (especially with GMG) and I want to be able to test Bicgstab for those models when I encounter them again to see if we accidentally have a non-symmetric operation in GMG somewhere (@tjhei and @jdannberg agreed this is a good idea).

[tjhei]

BiCGStab also works if the system is indefinite and not positive definite (in contrast to CG) as far as I know. That might be another reason to make this switch. Is this worth mentioning?

[gassmoeller]

/rebuild

[gassmoeller]

I think I addressed all comments. I put the check for asymmetric A block into a helper function to avoid the duplication and have added a note to say that Bicgstab can also handle indefinite matrices. Is there anything else left to do for this PR?

Just a quick reference for the future: In my tests GMRES was actually somewhat faster than Bicgstab for all the test cases I could find. Bicgstab reduces the number of iterations compared to CG, but it also needs 2 matrix-vector products per iteration instead of 1 (which is what gmres and cg need). GMRES would become slow if it needs many iterations because of the many scalar products to minimize the residual, but for the A block preconditioner the tolerance is low, and so GMRES rarely needs more than 10 iterations or so, and 10 scalar products seems to be cheaper than one matrix-vector product. However, this may change for larger models, or for harder problems, so let's go with Bicgstab.