SNES pseudo-timestepping fixes & improvements

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_<var name>``. 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                            |

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() {
@@ -1025,15 +1039,60 @@ 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 >= pseudo_squash_failure_threshold) {
+            // Squash variation in timestep between cells.
+            // pseudo_squash_lambda determines how much variation to keep.
+
+            switch (pseudo_squash_method) {
+            case BoutPseudoSquashMethod::affine:
+              // Modify dt_vec using Affine squash
+              // 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 = 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;
+              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, pseudo_squash_lambda);
+              }
+              PetscCall(VecRestoreArray(dt_vec, &dt_data));
+
+              pseudo_timestep =
+                  timestep * pow(pseudo_timestep / timestep, pseudo_squash_lambda);
+              break;
+            };
+
+            // Anti-windup: Calculate the effective alpha parameter
+            pseudo_alpha = mean(local_residual * pseudo_timestep, 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
@@ -1502,7 +1561,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()]);
         }
 

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,