From d9db28291f941d9ec74ec84b7dc7d246c6852d7c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 16 Aug 2022 15:11:43 +0100 Subject: [PATCH 1/3] Use Options& in LaplaceNaulin and document the options --- src/invert/laplace/impls/naulin/naulin_laplace.cxx | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index c722661b14..6d5731484d 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -153,6 +153,7 @@ LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc, Mesh *mesh_in) delp2solver(nullptr), naulinsolver_mean_its(0.), ncalls(0) { ASSERT1(opt != nullptr); // An Options pointer should always be passed in by LaplaceFactory + Options& options = *opt; Acoef.setLocation(location); C1coef.setLocation(location); @@ -160,10 +161,12 @@ LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc, Mesh *mesh_in) Dcoef.setLocation(location); // Get options - OPTION(opt, rtol, 1.e-7); - OPTION(opt, atol, 1.e-20); - OPTION(opt, maxits, 100); - OPTION(opt, initial_underrelax_factor, 1.); + rtol = options["rtol"].doc("relative tolerance").withDefault(1.e-7); + atol = options["atol"].doc("absolute tolerance").withDefault(1.e-20); + maxits = options["maxits"].doc("maximum number of iterations").withDefault(100); + initial_underrelax_factor = options["initial_underrelax_factor"] + .doc("Initial underrelaxation factor for the fixed point iteration.") + .withDefault(1.0); ASSERT0(initial_underrelax_factor > 0. and initial_underrelax_factor <= 1.); delp2solver = create(opt->getSection("delp2solver"), location, localmesh); std::string delp2type; From abaf1f96c0cde6807666cc75146d63bfaaad19ee Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 16 Aug 2022 15:28:48 +0100 Subject: [PATCH 2/3] Make iteration more robust and give more options in LaplaceNaulin Previously a cycle decreasing the underrelax_factor was triggered if the error increased on an iteration. However sometimes it can be more robust to just continue the iteration with moderate increases in the error - there are cases when decreasing the underrelax_factor never results in a decreasing error, but continuing the main iteration loop for a few steps would. This commit tries to handle this kind of case, with several new parameters which can be adjusted to control how and when underrelax_factor is decreased. --- .../laplace/impls/naulin/naulin_laplace.cxx | 37 ++++++++++++++++++- .../laplace/impls/naulin/naulin_laplace.hxx | 17 +++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index 6d5731484d..9e4634e748 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -168,6 +168,25 @@ LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc, Mesh *mesh_in) .doc("Initial underrelaxation factor for the fixed point iteration.") .withDefault(1.0); ASSERT0(initial_underrelax_factor > 0. and initial_underrelax_factor <= 1.); + underrelax_threshold = options["underrelax_threshold"] + .doc("Threshold for relative increase of error in a step that triggers decrease of " + "the underrelaxation factor.") + .withDefault(1.5); + ASSERT0(underrelax_threshold >= 1.); + underrelax_decrease_factor = options["underrelax_decrease_factor"] + .doc("Factor to decrease underrelax_factor at each stage of the sub-loop if " + "underrelax_threshold was crossed.") + .withDefault(0.9); + ASSERT0(underrelax_decrease_factor < 1. && underrelax_decrease_factor > 0.); + underrelax_decrease_maxits = options["underrelax_decrease_maxits"] + .doc("Maximum number of iterations in the decreasing-underrelax_factor subcycle " + "before trying to continue the main iteration loop.") + .withDefault(10); + underrelax_recovery = options["underrelax_recovery"] + .doc("Factor to increase underrelax_factor by at the end of a successful iteration " + "if it has been decreased below initial_underrelax_factor.") + .withDefault(1.1); + ASSERT0(underrelax_recovery >= 1.); delp2solver = create(opt->getSection("delp2solver"), location, localmesh); std::string delp2type; opt->getSection("delp2solver")->get("type", delp2type, "cyclic"); @@ -297,9 +316,10 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { throw BoutException("LaplaceNaulin error: Not converged within maxits=%i iterations.", maxits); } - while (error_abs > last_error) { + int local_count = 0; + while (error_abs > last_error * underrelax_threshold) { // Iteration seems to be diverging... try underrelaxing and restart - underrelax_factor *= .9; + underrelax_factor *= underrelax_decrease_factor; ++underrelax_count; // Restart from b_x_pair_old - that was our best guess @@ -312,6 +332,12 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { error_abs = max(abs(error3D, "RGN_NOBNDRY"), true, "RGN_NOBNDRY"); error_rel = error_abs / RMS_rhsOverD; + ++local_count; + if (local_count > underrelax_decrease_maxits) { + // Give up on trying to underrelax. Attempt to continue iteration anyway... + break; + } + // effectively another iteration, so increment the counter ++count; if (count>maxits) { @@ -322,6 +348,13 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { // Might have met convergence criterion while in underrelaxation loop if (error_rel initial_underrelax_factor) { + underrelax_factor = initial_underrelax_factor; + } + } + last_error = error_abs; b_x_pair_old = b_x_pair; b_x_pair = calc_b_x_pair(underrelax_factor*bnew + (1. - underrelax_factor)*b_x_pair.first, b_x_pair.second); diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.hxx b/src/invert/laplace/impls/naulin/naulin_laplace.hxx index 103364b8be..5e6c610668 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.hxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.hxx @@ -140,6 +140,23 @@ private: /// less than or equal to 1. Value of 1 means no underrelaxation BoutReal initial_underrelax_factor{1.}; + /// Factor by which the error may increase since the previous iteration before + /// triggering a sub-cycle with decreasing underrelax_factor to try and prevent + /// divergence. + BoutReal underrelax_threshold{1.5}; + + /// Factor by which to decrease underrelax_factor at each stage of the sub-loop + /// triggered if underrelax_threshold is crossed. + BoutReal underrelax_decrease_factor{0.9}; + + /// Maximum number of iterations in the decreasing-underrelax_factor subcycle before + /// trying to continue the main iteration loop. + int underrelax_decrease_maxits{10}; + + /// If underrelax_factor has been decreased, increase it again by this factor at the end + /// of a successful iteration to try and speed up subsequent convergence. + BoutReal underrelax_recovery{1.1}; + /// Mean number of iterations taken by the solver BoutReal naulinsolver_mean_its; From d78e861084ee3db1f1ed3267f99ffe5be6e24ffb Mon Sep 17 00:00:00 2001 From: johnomotani <3958036+johnomotani@users.noreply.github.com> Date: Fri, 17 Jan 2025 10:03:35 +0000 Subject: [PATCH 3/3] Apply clang-format changes --- .../laplace/impls/naulin/naulin_laplace.cxx | 43 +++++++++++-------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index c67e06b970..e6f68d850d 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -171,28 +171,35 @@ LaplaceNaulin::LaplaceNaulin(Options* opt, const CELL_LOC loc, Mesh* mesh_in, options["atol_accept"].doc("Accept this atol after maxits").withDefault(atol); maxits = options["maxits"].doc("maximum number of iterations").withDefault(100); - initial_underrelax_factor = options["initial_underrelax_factor"] - .doc("Initial underrelaxation factor for the fixed point iteration.") - .withDefault(1.0); + initial_underrelax_factor = + options["initial_underrelax_factor"] + .doc("Initial underrelaxation factor for the fixed point iteration.") + .withDefault(1.0); ASSERT0(initial_underrelax_factor > 0. and initial_underrelax_factor <= 1.); underrelax_threshold = options["underrelax_threshold"] - .doc("Threshold for relative increase of error in a step that triggers decrease of " - "the underrelaxation factor.") - .withDefault(1.5); + .doc("Threshold for relative increase of error in a step " + "that triggers decrease of " + "the underrelaxation factor.") + .withDefault(1.5); ASSERT0(underrelax_threshold >= 1.); - underrelax_decrease_factor = options["underrelax_decrease_factor"] - .doc("Factor to decrease underrelax_factor at each stage of the sub-loop if " - "underrelax_threshold was crossed.") - .withDefault(0.9); + underrelax_decrease_factor = + options["underrelax_decrease_factor"] + .doc("Factor to decrease underrelax_factor at each stage of the sub-loop if " + "underrelax_threshold was crossed.") + .withDefault(0.9); ASSERT0(underrelax_decrease_factor < 1. && underrelax_decrease_factor > 0.); - underrelax_decrease_maxits = options["underrelax_decrease_maxits"] - .doc("Maximum number of iterations in the decreasing-underrelax_factor subcycle " - "before trying to continue the main iteration loop.") - .withDefault(10); - underrelax_recovery = options["underrelax_recovery"] - .doc("Factor to increase underrelax_factor by at the end of a successful iteration " - "if it has been decreased below initial_underrelax_factor.") - .withDefault(1.1); + underrelax_decrease_maxits = + options["underrelax_decrease_maxits"] + .doc( + "Maximum number of iterations in the decreasing-underrelax_factor subcycle " + "before trying to continue the main iteration loop.") + .withDefault(10); + underrelax_recovery = + options["underrelax_recovery"] + .doc("Factor to increase underrelax_factor by at the end of a successful " + "iteration " + "if it has been decreased below initial_underrelax_factor.") + .withDefault(1.1); ASSERT0(underrelax_recovery >= 1.); delp2solver = create(opt->getSection("delp2solver"), location, localmesh); std::string delp2type;