defect correction + AMG convergence issue

[tjhei]

As reported by @naliboff : Convergence behavior of a model with AMG an defect correction takes many more iterations to converge on current main compared to 2.5. Let's investigate. continental_extension_iterated-defect-correction-stokes.prm.txt

[tjhei]

I can reproduce and I am running git bisect...

[tjhei]

7955a55b25f9b32fdbd83930756cceeef60608e7 is the first bad commit
commit 7955a55b25f9b32fdbd83930756cceeef60608e7
Author: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date:   Fri Jul 14 16:25:40 2023 -0700

    Be safer and more correct when using strain rate

 source/material_model/rheology/visco_plastic.cc    |  18 +-
 source/material_model/visco_plastic.cc             |   6 +-
 source/simulator/assembly.cc                       |   1 +
 source/simulator/melt.cc                           |   1 +
 tests/continental_extension/screen-output          |   4 +-
 tests/visco_plastic_derivatives_2d.prm             |   1 +
 tests/visco_plastic_derivatives_2d/screen-output   | 112 +++++------
 tests/visco_plastic_derivatives_3d.prm             |   1 +
 tests/visco_plastic_derivatives_3d/screen-output   | 213 ++++++++++-----------
 .../screen-output                                  |   4 +-
 10 files changed, 186 insertions(+), 175 deletions(-)

[tjhei]

see #5239

[bangerth]

@gassmoeller That's one of yours.

[naliboff]

@bangerth @tjhei @gassmoeller - I looked over #5329 again, and I still think all of the changes in there are fundamentally correct.

What it is not clear to me is why these changes are having such a large affect on the number of required defection correction Stokes iterations, as reported and discussed on the user forums, while the number of iterations for the iterated Stokes solver scheme does not change at all.

@MFraters @bangerth - Are there cases where it conceptually makes sense that the regular picard iterations would converge much more quickly than defect correction picard iterations? If yes and the changes to the tests in #5239 still seem correct, perhaps the last check would be to rerun the full newton solver benchmark suite?

[anne-glerum]

@naliboff Looking at #5239, I do think it introduced a change in the used strain rate in some cases that is not in line with the user-set minimum and reference strain rate values.

Before:

 const bool use_reference_strainrate = (this->get_timestep_number() == 0) &&
                                       (in.strain_rate[i].norm() <= std::numeric_limits<double>::min());

After:

const bool use_reference_strainrate = this->simulator_is_past_initialization() == false
                                       ||
                                       (this->get_timestep_number() == 0 &&
                                        this->get_nonlinear_iteration() == 0)
                                       ||
                                       (in.strain_rate[i].norm() <= std::numeric_limits<double>::min());

In the 'before' case, the reference_strainrate was effectively only used during the first timestep and probably only during the first nonlinear iteration. Note that there is a second parameter min_strain_rate that takes care of strain rates being too low later on in the code. The reference strain rate was only intended to provide a user-set value for the first non-linear iteration based on for example the prescribed velocity boundary conditions.

In the 'after' case, the reference strain rate will also be used in case the strain rate is very low (ie the case that was previously caught by the min_strain_rate value). By default this ref_strain_rate is 5 orders of magnitude larger than min_strain_rate.

I don't know whether this change would lead to bad convergence, but I do think it should be reverted.

[naliboff]

@anne-glerum - Thanks for catching this and per our discussion, we were all in agreement that this should be reverted.