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

[YiminJin]

@naliboff I agree with the speculation:

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

Here is a brief explanation:

The figure above shows the return mapping procedure on $\pi$-plane. The physical meanings of the vectors are:
(1) $\vec{OA} = 2\eta\varepsilon = \tau^{trial}$,
(2) $\vec{AB} = 2\eta\delta\varepsilon = \delta\tau^{trial}$,
(3) $\vec{AC} \approx 2\eta(\delta\varepsilon - n\otimes n:\delta\varepsilon)$,
(4) $\vec{ED} \approx 2\eta^{\rm eff}(\delta\varepsilon - n\otimes n:\delta\varepsilon) = \delta\tau$,

where $n := \varepsilon / \Vert\varepsilon\Vert$ = $\vec{OA}/\Vert\vec{OA}\Vert$. It should be noted that in (3) and (4) we have made the approximation $\vec{OB}/\Vert\vec{OB}\Vert \approx n$, and the approximation is reasonable only when $\delta\varepsilon$ is small (as shown in Fig. (a)).

However, at the first time step the initial guess is greatly deviated from the true solution, and $\delta\varepsilon$ is large (remember that $\delta\varepsilon$ is not essentially a small value, it is determined by r/r'). In this case, (3) and (4) are bad approximations, and we can see that the resulting stress is far from the yield surface (which means that the consistent condition $\delta F = 0$ is violated, as shown in Fig. (b)). Hence, we need to scale the Newton increment in order to keep the resulting stress close to the yielding surface (point D'), and that is why we need the Newton residual scaling method to ensure convergence at the first time step.

I think this also explains why the Newton solver performs better with VEP models: when elasticity is included, plastic yielding does not occur at the first time step in most cases (or only occurs at a few points), so $\delta\varepsilon$ is small.

[YiminJin]

Leaving aside the previous comment, I recommend an alternative method to initialize the plastic strain on particles. Here is my plugin:
initial_plastic_strain.cc.txt
This plugin reinitializes particle properties with the signal slot post_set_initial_state. I generate a random value for each cell, which is applied to all the particles inside the cell. Then if we change the interpolation scheme to "cell average", the results look like (GMG, Iterated Advection and Newton Stokes, $10^{-7}$):

However, in order to make this plugin work, I have to modify several lines in ASPECT: the function Particle::World<dim>::setup_initial_state() is connected to the signal slot post_set_initial_state itself, which will overwrite the function in my plugin. So I need to remove it from Particle::World<dim>::connect_to_signals() and call it explicitly in Simulator<dim>::run(). I do not know if there are special reasons to treat some of the important functions in Particle::World<dim> as signals. If not, I suggest calling them explicitly in the simulator, which would provide more freedom to the users. @gassmoeller

[gassmoeller]

So I need to remove it from Particle::World::connect_to_signals() and call it explicitly in Simulator::run(). I do not know if there are special reasons to treat some of the important functions in Particle::World as signals.

I do not see a reason why the function has to be connected to the signal. I think we implemented this at a time when the particle handler was rarely used, and only seen as one of the possible postprocessors, so we didnt think it to be appropriate to put the call into the main simulator class. This has changed since then, and particle worlds are now owned by the simulator directly. So in short, I would support removing the use of the signal, and instead calling Particle::World<dim>::setup_initial_state() directly from the place before the signal is triggered.

Another option would be of course to modify the particle plugin initialization to allow your routine as one of the options (computing a random value per cell and then applying that to the particles as initial strain). Then you wouldnt need to make any change outside the particle plugin.

[naliboff]

@YiminJin - Thanks for the explanation above regarding the Newton method and divergence that can potentially arise from the first guess of the strain rate. Indeed it makes much more sense now why the Newton solver performs much better with VEP models.

As you noted, the last series of models is so sensitive to the initial strain rate due the element-scale randomization of the initial plastic strain. Let's discuss further with @MFraters and others to sort out a path forward. Thank you again for laying out the helpful diagrams and equations!