From 9c262818009ada6eab266654a899da34a6120b71 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 11 Sep 2026 13:46:42 -0700 Subject: [PATCH 1/4] SNESSolver::updatePseudoTimestep fix typo Copy-paste error in calculation of neighboring timesteps. --- src/solver/impls/snes/snes.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index eb786822a3..403f9217b9 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1502,7 +1502,7 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() { if (i3d.y() != 0) { min_neighboring_dt = std::min(min_neighboring_dt, pseudo_timestep[i3d.ym()]); } - if (i3d.x() != mesh->LocalNy - 1) { + if (i3d.y() != mesh->LocalNy - 1) { min_neighboring_dt = std::min(min_neighboring_dt, pseudo_timestep[i3d.yp()]); } From fd8dc16ea8c4d923f380427f7086de545a7444ea Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 11 Sep 2026 13:48:56 -0700 Subject: [PATCH 2/4] SNESSolver: Improve failure handling in pseudo_transient mode Modify dt_vec PETSc Vec, pseudo_timestep Field3D and pseudo_alpha on SNES failure. When dt_vec was modified but pseudo_alpha was not, the PID controller would become saturated so that timestep increases were determined by caps on increases. --- src/solver/impls/snes/snes.cxx | 46 +++++++++++++++++++++++++++++++--- 1 file changed, 43 insertions(+), 3 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 403f9217b9..f9c0616187 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1025,15 +1025,55 @@ int SNESSolver::run() { if (snes_failures == max_snes_failures - 1) { // Last chance. Set to uniform smallest timestep PetscCall(VecSet(dt_vec, dt_min_reset)); + pseudo_timestep = dt_min_reset; + + // Scale down pseudo_alpha so that PID controller isn't saturated + pseudo_alpha = pseudo_alpha_minimum; + + } else if (snes_failures >= 5) { + // Squash variation in timestep + // dt_vec <- lambda * dt_vec + (1 - lambda) * timestep + const BoutReal lambda = 0.5; + const bool affine_squash = false; + + if (affine_squash) { + // Modify dt_vec using Affine squash + PetscCall(VecScale(dt_vec, lambda)); + PetscCall(VecShift(dt_vec, (1.0 - lambda) * timestep)); + + // Modify pseudo_timestep + pseudo_timestep = lambda * pseudo_timestep + (1 - lambda) * timestep; + } else { + // Log squash + // dt_vec <- timestep * (dt_vec / timestep)^lambda + PetscInt size; + PetscCall(VecGetLocalSize(dt_vec, &size)); + BoutReal* dt_data = nullptr; + PetscCall(VecGetArray(dt_vec, &dt_data)); + for (PetscInt i = 0; i != size; ++i) { + dt_data[i] = timestep * std::pow(dt_data[i] / timestep, lambda); + } + PetscCall(VecRestoreArray(dt_vec, &dt_data)); + + pseudo_timestep = timestep * pow(pseudo_timestep / timestep, lambda); + } + + // Anti-windup: Calculate the effective alpha parameter + Field3D pseudo_alpha_effective = local_residual * pseudo_timestep; + pseudo_alpha = mean(pseudo_alpha_effective, true); - } else if (snes_failures == 5) { - // Set uniform timestep - PetscCall(VecSet(dt_vec, timestep)); } else { // Global scaling of timesteps // Note: A better strategy might be to reduce timesteps // in problematic cells. PetscCall(VecScale(dt_vec, timestep_factor_on_failure)); + + pseudo_timestep *= timestep_factor_on_failure; + + // Scale alpha down by the same amount + // If this is not done then PID controller can 'wind up' + // because the controller keeps increasing alpha between failures. + pseudo_alpha *= timestep_factor_on_failure; } } else { // Try a smaller timestep From ba9bdab09bf9894618d28048be3f30567c987ebe Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 18 Sep 2026 09:40:41 -0700 Subject: [PATCH 3/4] SNESSolver: Add psuedo_squash settings One cause of solver failure can be very large variation of the timestep between cells. When SNES repeatedly fails the pseudo timestep can therefore be "squashed" to reduce variation. This adds options to tune this process: - `pseudo_squash_failure_threshold` is an integer (default 5). Once the number of snes failures reaches this threshold squashing is performed (subsequent failures will be squashed again). When snes succeeds the counter is reset to zero. - `pseudo_squash_method` is an enum, `affine` or `log` to determine whether the squashing is linear or logarithmic. - `pseudo_squash_lambda` is a float between 0 and 1. It determines how much of the variation to keep. --- src/solver/impls/snes/snes.cxx | 51 +++++++++++++++++++++++----------- src/solver/impls/snes/snes.hxx | 11 +++++++- 2 files changed, 45 insertions(+), 17 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index f9c0616187..c4242ec5d9 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -104,7 +104,7 @@ PetscErrorCode withOptionalSubvectors(Func operation, IS indices, Args... args) } for (std::size_t i = acquired; i > 0; --i) { - PetscErrorCode restore_ierr = + const PetscErrorCode restore_ierr = VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]); if (ierr == PETSC_SUCCESS) { ierr = restore_ierr; @@ -332,6 +332,18 @@ SNESSolver::SNESSolver(Options* opts) pseudo_max_ratio((*options)["pseudo_max_ratio"] .doc("PTC maximum timestep ratio between neighbors") .withDefault(2.)), + pseudo_squash_failure_threshold( + (*options)["pseudo_squash_failure_threshold"] + .doc("Squash timestep variation when snes failures reaches this threshold") + .withDefault(5)), + pseudo_squash_method( + (*options)["pseudo_squash_method"] + .doc("Method to apply when squashing pseudo timesteps: affine or log.") + .withDefault(BoutPseudoSquashMethod::log)), + pseudo_squash_lambda((*options)["pseudo_squash_lambda"] + .doc("How much variation to keep? 0 = No variation; 1 = " + "Full variation (no squashing)") + .withDefault(0.5)), timestep_control((*options)["timestep_control"] .doc("Timestep control method") .withDefault(BoutSnesTimestep::pid_nonlinear_its)), @@ -427,6 +439,8 @@ SNESSolver::SNESSolver(Options* opts) .doc("Which Jacobian to save: system, scaled, or rhs") .withDefault(bout::JacobianExportKind::system)) { supports_constraints = true; // This solver can handle constraints + + ASSERT0((pseudo_squash_lambda >= 0.0) and (pseudo_squash_lambda <= 1.0)); } SNESSolver::~SNESSolver() { @@ -1030,20 +1044,23 @@ int SNESSolver::run() { // Scale down pseudo_alpha so that PID controller isn't saturated pseudo_alpha = pseudo_alpha_minimum; - } else if (snes_failures >= 5) { - // Squash variation in timestep - // dt_vec <- lambda * dt_vec + (1 - lambda) * timestep - const BoutReal lambda = 0.5; - const bool affine_squash = false; + } else if (snes_failures >= pseudo_squash_failure_threshold) { + // Squash variation in timestep between cells. + // pseudo_squash_lambda determines how much variation to keep. - if (affine_squash) { + switch (pseudo_squash_method) { + case BoutPseudoSquashMethod::affine: // Modify dt_vec using Affine squash - PetscCall(VecScale(dt_vec, lambda)); - PetscCall(VecShift(dt_vec, (1.0 - lambda) * timestep)); + // dt_vec <- lambda * dt_vec + (1 - lambda) * timestep + + PetscCall(VecScale(dt_vec, pseudo_squash_lambda)); + PetscCall(VecShift(dt_vec, (1.0 - pseudo_squash_lambda) * timestep)); // Modify pseudo_timestep - pseudo_timestep = lambda * pseudo_timestep + (1 - lambda) * timestep; - } else { + pseudo_timestep = pseudo_squash_lambda * pseudo_timestep + + (1 - pseudo_squash_lambda) * timestep; + break; + case BoutPseudoSquashMethod::log: // Log squash // dt_vec <- timestep * (dt_vec / timestep)^lambda PetscInt size; @@ -1051,16 +1068,18 @@ int SNESSolver::run() { BoutReal* dt_data = nullptr; PetscCall(VecGetArray(dt_vec, &dt_data)); for (PetscInt i = 0; i != size; ++i) { - dt_data[i] = timestep * std::pow(dt_data[i] / timestep, lambda); + dt_data[i] = + timestep * std::pow(dt_data[i] / timestep, pseudo_squash_lambda); } PetscCall(VecRestoreArray(dt_vec, &dt_data)); - pseudo_timestep = timestep * pow(pseudo_timestep / timestep, lambda); - } + pseudo_timestep = + timestep * pow(pseudo_timestep / timestep, pseudo_squash_lambda); + break; + }; // Anti-windup: Calculate the effective alpha parameter - Field3D pseudo_alpha_effective = local_residual * pseudo_timestep; - pseudo_alpha = mean(pseudo_alpha_effective, true); + pseudo_alpha = mean(local_residual * pseudo_timestep, true); } else { // Global scaling of timesteps diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 57a4c4e842..40d39cf059 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -74,6 +74,10 @@ BOUT_ENUM_CLASS(BoutSnesOutput, fixed_time_interval, ///< Output at fixed time intervals residual_ratio); ///< When the residual is reduced by a given ratio +BOUT_ENUM_CLASS(BoutPseudoSquashMethod, + affine, ///< Affine dt_vec <- lambda * dt_vec + (1 - lambda) * timestep + log); ///< Log squash dt_vec <- timestep * (dt_vec / timestep)^lambda + /// Uses PETSc's SNES interface to find a steady state solution to a /// nonlinear ODE by integrating in time with Backward Euler class SNESSolver : public Solver { @@ -189,7 +193,12 @@ private: 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 - Vec dt_vec; ///< Each quantity can have its own timestep + int pseudo_squash_failure_threshold; ///< Squash timestep variation when snes failures exceed this threshold + BoutPseudoSquashMethod + pseudo_squash_method; ///< Method to apply when squashing pseudo timesteps + BoutReal + pseudo_squash_lambda; ///< How much variation to keep? 0 = No variation; 1 = Full variation (no squashing). + Vec dt_vec; ///< Each quantity can have its own timestep /// Adjust the global timestep BoutReal updateGlobalTimestep(BoutReal timestep, int nl_its, From 2b5183b46915e87e241e9c03426d3bf0b274ca9f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 18 Sep 2026 13:09:19 -0700 Subject: [PATCH 4/4] SNESSolver: Update PTC method docs Document the new timestep squashing options. Simplify descriptions of the method and timestep adjustment in the PTC method. --- manual/sphinx/user_docs/time_integration.rst | 129 ++++++++++++++----- 1 file changed, 94 insertions(+), 35 deletions(-) diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index ff3d5942ae..b72f47bcd1 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -706,9 +706,9 @@ Pseudo-Transient Continuation and Switched Evolution Relaxation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When ``equation_form = pseudo_transient`` the solver uses -Pseudo-Transient Continuation (PTC). This is a robust numerical -technique for solving steady-state problems that are too nonlinear for -direct Newton iteration. Instead of solving the steady-state system +Pseudo-Transient Continuation (PTC). This method helps with steady +state problems that are too nonlinear for direct Newton iteration. +Instead of solving the steady-state system **F(u) = 0** directly, PTC solves a modified time-dependent problem: .. math:: @@ -719,27 +719,32 @@ where :math:`\tau` is a pseudo-time variable (not physical time) and :math:`M(u) is a preconditioning matrix. As :math:`\tau \to \infty`, the solution converges to the steady state **F(u) = 0**. -The key advantage of PTC is that it transforms a difficult root-finding problem -into a sequence of easier initial value problems. Poor initial guesses that would -cause Newton's method to diverge can still reach the solution via a stable -pseudo-transient path. +PTC turns one hard root-finding problem into a set of easier time +steps. A poor first guess can still reach the steady state. The Switched Evolution Relaxation (SER) method is a spatially adaptive -variant of PTC that allows each cell to use a different -pseudo-timestep :math:`\Delta\tau_i`. The timestep in each cell adapts -based on the local residual, allowing the algorithm to take large -timesteps in well-behaved regions (fast convergence), while taking -small timesteps in difficult regions (stable advancement). The the -same :math:`\Delta\tau_i` is used for all equations (density, -momentum, energy etc.) within each cell. This maintains coupling -between temperature, pressure, and composition through the equation of -state. +form of PTC. Each cell can use its own pseudo-timestep +:math:`\Delta\tau_i`. The timestep in each cell changes with the local +residual. Cells that behave well can take large steps, while cells that are +hard to solve take small steps. The same :math:`\Delta\tau_i` is used +for all equations in one cell, maintaining the equation of state within +each cell. **Key parameters:** ``pseudo_max_ratio`` (default: 2.0) - Maximum allowed ratio of timesteps between neighboring cells. This prevents - sharp spatial gradients in convergence rate. + Largest allowed ratio of timesteps between nearby cells. + +``pseudo_squash_failure_threshold`` (default: 5) + Start to squash the spread in local pseudo-timesteps after this many + SNES failures in a row. + +``pseudo_squash_method`` (default: ``log``) + How to squash the local pseudo-timesteps. Use ``affine`` or ``log``. + +``pseudo_squash_lambda`` (default: 0.5) + How much of the old spread to keep. ``0`` gives one common timestep. + ``1`` keeps the old spread. **Example PTC configuration:** @@ -754,6 +759,9 @@ state. # SER parameters timestep_control = pid_nonlinear_its # Scale timesteps based on iterations pseudo_max_ratio = 2.0 # Limit neighbor timestep ratio + pseudo_squash_failure_threshold = 5 + pseudo_squash_method = log + pseudo_squash_lambda = 0.5 # Tolerances atol = 1e-7 @@ -778,17 +786,17 @@ is computed as: \Delta\tau_i = \frac{\alpha}{||R_i||} -Larger values allow more aggressive timestepping. The default is to use -a fixed ``pseudo_alpha`` but a better strategy is to enable the PID controller -that adjusts this parameter based on the nonlinear solver convergence. +Large values give larger timesteps. By default ``pseudo_alpha`` is +fixed, but you can also let the PID controller change it based on the +nonlinear solve history. The timestep is limited to be between ``dt_min_reset`` and ``max_timestep``. In addition the timestep is limited between 0.67 × previous timestep and 1.5 × previous timestep, to limit sudden changes in timestep. -In practice this strategy seems to work well, though problems could -arise when residuals become very small. +This often works well, but very small residuals can still cause +problems. **history_based** @@ -816,8 +824,8 @@ become small the method switches to ``history_based``. PID Controller ^^^^^^^^^^^^^^ -When using the PTC method the PID controller can be used to dynamically -adjust ``pseudo_alpha`` depending on the nonlinearity of the system: +When you use PTC, the PID controller can change ``pseudo_alpha`` to +adjust to the nonlinearity of the system: .. code-block:: ini @@ -828,15 +836,43 @@ adjust ``pseudo_alpha`` depending on the nonlinearity of the system: kI = 0.3 # Integral gain kD = 0.2 # Derivative gain -The PID controller adjusts ``pseudo_alpha``, scaling all cell -timesteps together, to maintain approximately ``target_its`` nonlinear +The PID controller adjusts ``pseudo_alpha``. This scales all cell +timesteps together and aims for about ``target_its`` nonlinear iterations per solve. -With this enabled the solver uses the number of nonlinear iterations -to scale timesteps globally, and residuals to scale timesteps locally. +With this on, the solver uses the number of nonlinear iterations to +scale timesteps for the whole domain, and uses residuals to scale +timesteps in each cell. + +On repeated SNES failures, the solver now also scales down or resets +``pseudo_alpha`` so the PID controller does not keep pushing the +timestep back up. + +Current limit: this anti-windup step uses only the ``Field3D`` local +residual and pseudo-timestep data. It does not include ``Field2D`` +parts. In mixed ``Field2D``/``Field3D`` cases, the reset reflects only +the ``Field3D`` part. In ``Field2D``-only cases, it does not give a +useful value. + Note that the PID controller has no effect on the ``history_based`` strategy because that strategy does not use ``pseudo_alpha``. +PTC failure handling +^^^^^^^^^^^^^^^^^^^^ + +If SNES fails in ``pseudo_transient`` mode, the solver first scales all +local pseudo-timesteps down by ``timestep_factor_on_failure``. + +After ``pseudo_squash_failure_threshold`` failures in a row, the solver +also squashes the spread in local pseudo-timesteps. This pulls them +toward the current global timestep. Use ``pseudo_squash_method`` to +pick the squash rule and ``pseudo_squash_lambda`` to set how strong the +squash is. + +On the last retry before the solver stops, it sets all local +pseudo-timesteps to ``dt_min_reset`` and sets ``pseudo_alpha`` to +``pseudo_alpha_minimum``. + Jacobian Finite Difference with Coloring ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1012,14 +1048,18 @@ Diagnostics and Monitoring diagnose = true # Print iteration info to screen diagnose_failures = true # Detailed diagnostics on failures -When ``equation_form = pseudo_transient``, the solver saves additional diagnostic fields: +When ``equation_form = pseudo_transient``, the solver saves extra diagnostic fields: -- ``snes_pseudo_residual``: Local residual in each cell +- ``snes_local_residual``: Local residual in each cell +- ``snes_global_residual``: Global RMS residual - ``snes_pseudo_timestep``: Local pseudo-timestep in each cell - ``snes_pseudo_alpha``: Global timestep scaling -These can be visualized to understand convergence behavior and identify -problematic regions. +These can help you see why the solve is slow or where it fails. + +The anti-windup update for ``snes_pseudo_alpha`` has one limit at +present: it uses only ``Field3D`` local residual and pseudo-timestep +data. It does not include ``Field2D`` parts. The residuals from the last nonlinear solve are also saved with names ``resid_``. Plotting these can help to understand which @@ -1031,11 +1071,30 @@ Summary of solver options +---------------------------+---------------+----------------------------------------------------+ | Option | Default |Description | +===========================+===============+====================================================+ -| pseudo_time | false | Pseudo-Transient Continuation (PTC) method, using | -| | | a different timestep for each cell. | +| equation_form | rearranged_ | Choose the SNES solve form. Use | +| | backward_ | ``pseudo_transient`` for PTC. | +| | euler | | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_alpha | 100*atol*dt | Sets local timestep in ``inverse_residual`` mode | +| | | with ``dt = pseudo_alpha / residual`` | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_alpha_minimum | 0.1*pseudo_ | Smallest allowed value for ``pseudo_alpha`` | +| | alpha | | +---------------------------+---------------+----------------------------------------------------+ | pseudo_max_ratio | 2. | Maximum timestep ratio between neighboring cells | +---------------------------+---------------+----------------------------------------------------+ +| pseudo_growth_factor | 1.1 | Growth factor in ``history_based`` mode | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_reduction_factor | 0.5 | Reduction factor in ``history_based`` mode | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_squash_failure_ | 5 | Start to squash local pseudo-timesteps after this | +| threshold | | many SNES failures in a row | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_squash_method | log | How to squash local pseudo-timesteps: ``affine`` | +| | | or ``log`` | ++---------------------------+---------------+----------------------------------------------------+ +| pseudo_squash_lambda | 0.5 | How much of the old timestep spread to keep | ++---------------------------+---------------+----------------------------------------------------+ | snes_type | newtonls | PETSc SNES nonlinear solver (try anderson, qn) | +---------------------------+---------------+----------------------------------------------------+ | ksp_type | gmres | PETSc KSP linear solver |