Solution divergence between different combinations of linear and nonlinear solver schemes

[naliboff]

This post is a follow-up from a related series of issues and pull requests (#4370, #4563, #6002) that discuss solver failure or solution divergence for different combinations linear (AMG, GMG) and nonlinear (iterated Stokes, iterated Newton Stokes, etc) solvers.

@anne-glerum, @egpuckett, and previously discussed these results a few weeks ago.

As a starting point, I first ran a modified version of benchmarks/viscoelastic_plastic_shear_bands/kaus_2010/kaus_2010_extension.prm with elasticity removed and other updates to the solver parameters based on recent testing. One examples of this PRM is attached.

This model was run for 0.01 Myr (100 time steps) was run for the following four parameter combinations:

1. Nonlinear solver tolerance = 1e-5, zero traction surface (no mesh deformation)
2. Nonlinear solver tolerance = 1e-7, zero traction surface (no mesh deformation)
3. Nonlinear solver tolerance = 1e-5, free surface (mesh deformation)
4. Nonlinear solver tolerance = 1e-7, free surface (mesh deformation)

For each of the four parameter combinations above, the following four solver combinations were tested:

1. Stokes solver type = block AMG, Nonlinear solver scheme = single Advection, iterated Stokes
2. Stokes solver type= block GMG, Nonlinear solver scheme = single Advection, iterated Stokes
3. Stokes solver type= block AMG, Nonlinear solver scheme = single Advection and iterated Newton Stokes
4. Stokes solver type= block GMG, Nonlinear solver scheme = single Advection and iterated Newton Stokes

At time = 0, the solutions (strain rate) are qualitatively very similar for a nonlinear solver tolerance of 1e-5 and zero traction surface or free surface:

Zero Traction Surface, t=0, Nonlinear Solver Tolerance = 1e-5

Free Surface, t=0 Myr, Nonlinear Solver Tolerance = 1e-5*

From t=0 to t=0.1 Myr, the two models (block AMG versus block GMG) using single Advection, Iterated Stokes still have qualitatively similar solutions, while the models using single Advection, iterated Newton Stokes diverge:

Zero Traction Surface, t=0.1 Myr, Nonlinear Solver Tolerance = 1e-5

Free Surface, t=0.1 Myr, Nonlinear Solver Tolerance = 1e-5*

Making the nonlinear solver tolerance stricter (1e-7) leads to qualitatively similar strain rate fields for both the zero traction and free surface cases:

Zero Traction Surface, t=0.1 Myr, Nonlinear Solver Tolerance = 1e-7

Free Surface, t=0.1 Myr, Nonlinear Solver Tolerance = 1e-7*

Notably, the results for the case of Nonlinear solver scheme = single Advection, iterated Stokes and a Nonlinear solver tolerance = 1e-5 or 1e-7 show little divergence.

Initial conclusions: For the same linear solver tolerance, a stricter nonlinear solver tolerance is needed for iterated Newton Stokes (and presumably iterated defect correction Stokes) related to iterated Stokes.

Next Steps:

1. Test to see if we can achieve qualitatively similar solutions for the continental extension cookbook with stricter solver tolerances across the four solver scheme combinations tested here.
2. Update the continental extension cookbook and related documentation accordingly.

@MFraters @bangerth @YiminJin @tjhei - This may be of interest.

Example PRM File:
kaus_2010_extension_gmg_iterated-newton-stokes.prm.txt

[naliboff]

An update from the post above with additional tests using a modified version of the continental extension cookbook.

Oddly, the solutions appear very similar at a nonlinear solver tolerance of 1e-5, but then diverge between nonlinear solver schemes for a tolerance of 1e-6. I'm not sure what to make of this.

The models would not converge to a nonlinear solver tolerance greater than 1e-6 even when the mix/max viscosity range was reduced or damper viscosity increased.

However, allowing the mesh on the left and right boundaries to deform (via additional tangential mesh boundary indicators) is producing more stable convergence in an additional set of tests, and I will provide an update once those tests have finished.

Time = 10 Myr, Nonlinear Solver Tolerance = 1e-5

Time = 10 Myr, Nonlinear Solver Tolerance = 1e-6

modified_continental_extension_cookbook.prm.txt

[tjhei]

Initial conclusions: For the same linear solver tolerance, a stricter nonlinear solver tolerance is needed for iterated Newton Stokes (and presumably iterated defect correction Stokes) related to iterated Stokes.

I am not quite sure I follow. I would pick a nonlinear solver tolerance that is accurate enough for your problem (if all schemes converge to 1e-12 nonlinear tolerance, I expect all solutions to be close to each other). After you pick your solver scheme and solver tolerance, you are free to choose a linear solver tolerance that a) makes it converge and b) gives a fast time to solution.

[naliboff]

Update results for when set Additional tangential mesh velocity boundary indicators = left, right is added to the extension model above, allowing convergence to a nonlinear solver tolerance of 1e-7.

iterated Advection and Stokes is still producing effectively the same solution, while the Newton solver is still diverging between solver tolerances and setup for this problem.

The next step is to simplify the problem a bit (use a single seed versus random seeds), as the the brick benchmark did show convergence between the four cases.

@MFraters @anne-glerum - I wonder if the divergence is due to the Newton solver using the reference strain as the first guess during each time step, rather than the velocity field from the last time step?

Time = 10 Myr, Nonlinear Solver Tolerance = 1e-7