From 5ec260f7fab2fb4f620fec9fce3ee24d5957f9d3 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 22 May 2026 09:22:03 -0700 Subject: [PATCH] snes: pseudo_alpha_minimum setting This is to prevent the PTC method from getting stuck when the alpha parameter gets very small, taking many iterations to recover. --- src/solver/impls/snes/snes.cxx | 42 ++++++++++++++++++++++------------ src/solver/impls/snes/snes.hxx | 1 + 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 31877f8aac..5b28ff413c 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -4,18 +4,25 @@ #include "snes.hxx" +#include +#include #include #include +#include +#include #include #include #include #include +#include #include +#include #include #include #include #include +#include #include #include "petscerror.h" @@ -306,7 +313,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { } // 3D fields for (int z = mesh->zstart; z <= mesh->zend; z++) { - int ind = ROUND(index(x, y, z)); + const int ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { PetscInt row = ind + i; @@ -342,7 +349,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { // 3D fields on this cell for (int j = 0; j < n3d; j++) { const PetscInt col = ind2 + j; - PetscErrorCode ierr = + const PetscErrorCode ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); if (ierr != PETSC_SUCCESS) { @@ -536,6 +543,9 @@ SNESSolver::SNESSolver(Options* opts) pseudo_alpha((*options)["pseudo_alpha"] .doc("Sets timestep using dt = alpha / residual") .withDefault(100. * atol * timestep)), + pseudo_alpha_minimum( + (*options)["pseudo_alpha_minimum"].doc("Minimum value of pseudo_alpha") + .withDefault(0.1 * pseudo_alpha)), pseudo_growth_factor((*options)["pseudo_growth_factor"] .doc("PTC growth factor on success") .withDefault(1.1)), @@ -884,7 +894,7 @@ int SNESSolver::run() { int saved_jacobian_lag = 0; int loop_count = 0; - BoutReal start_global_residual = global_residual; + const BoutReal start_global_residual = global_residual; do { if ((output_trigger == BoutSnesOutput::fixed_time_interval && (simtime >= target)) || (output_trigger == BoutSnesOutput::residual_ratio @@ -1181,8 +1191,10 @@ int SNESSolver::run() { if (equation_form == BoutSnesEquationForm::pseudo_transient) { // Adjust pseudo_alpha to globally scale timesteps - pseudo_alpha = updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate, - max_timestep * atol * 100); + pseudo_alpha = std::max({ + updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate, + max_timestep * atol * 100), + pseudo_alpha_minimum}); // Adjust local timesteps PetscCall(updatePseudoTimestepping()); @@ -1216,7 +1228,7 @@ int SNESSolver::run() { // Note: If simtime = target then alpha = 0 // and output_x = snes_x - BoutReal alpha = (simtime - target) / dt; + const BoutReal alpha = (simtime - target) / dt; // output_x <- alpha * x0 + (1 - alpha) * output_x VecAXPBY(output_x, alpha, 1. - alpha, x0); @@ -1437,7 +1449,7 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() { } if (count > 0) { // Adjust timestep for these quantities - BoutReal new_timestep = updatePseudoTimestep( + const BoutReal new_timestep = updatePseudoTimestep( dt_data[idx], local_residual_2d_prev[i2d], local_residual_2d[i2d]); for (int i = 0; i != count; ++i) { dt_data[idx++] = new_timestep; @@ -1529,9 +1541,10 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() { /// rapid changes in timestep. BoutReal SNESSolver::updatePseudoTimestep_inverse_residual(BoutReal previous_timestep, BoutReal current_residual) { - return std::min( - {std::max({pseudo_alpha / current_residual, previous_timestep / 1.5, dt_min_reset}), - 1.5 * previous_timestep, max_timestep}); + return std::max({ + std::min({std::max({pseudo_alpha / current_residual, previous_timestep / 1.5}), + 1.5 * previous_timestep, max_timestep}), + dt_min_reset}); } // Strategy based on history of residuals @@ -1552,7 +1565,8 @@ BoutReal SNESSolver::updatePseudoTimestep_history_based(BoutReal previous_timest if (reduction_ratio < 0.8) { return std::min(0.5 * (pseudo_growth_factor + 1.) * previous_timestep, max_timestep); - } else if (reduction_ratio > 1.2) { + } + if (reduction_ratio > 1.2) { return std::max(0.5 * (pseudo_reduction_factor + 1) * previous_timestep, dt_min_reset); } @@ -1780,10 +1794,8 @@ PetscErrorCode SNESSolver::scaleJacobian(Mat Jac_new) { // Calculate a norm of this row of the Jacobian PetscScalar norm = 0.0; for (int col = 0; col < ncols; col++) { - PetscScalar absval = std::abs(vals[col]); - if (absval > norm) { - norm = absval; - } + const PetscScalar absval = std::abs(vals[col]); + norm = std::max({norm, absval}); // Can we identify small elements and remove them? // so we don't need to calculate them next time } diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 40412f83b7..fbcf3baddf 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -153,6 +153,7 @@ private: // These are used if equation_form = pseudo_transient BoutPTCStrategy pseudo_strategy; ///< Strategy to use when setting timesteps BoutReal pseudo_alpha; ///< dt = alpha / residual + BoutReal pseudo_alpha_minimum; ///< Minimum value of alpha BoutReal pseudo_growth_factor; ///< Timestep increase 1.1 - 1.2 BoutReal pseudo_reduction_factor; ///< Timestep decrease 0.5 BoutReal pseudo_max_ratio; ///< Maximum timestep ratio between neighboring cells