From dc3720f7b57aeb9cf0dde05374dffd664f0e37df Mon Sep 17 00:00:00 2001
From: Ben Dudson <dudson2@llnl.gov>
Date: Fri, 10 Oct 2025 11:36:51 -0700
Subject: [PATCH] snes solver: Consider timestep change at output

Previously the timestep adjustment was inside a conditional
`if (looping) {...`. The variable `looping` is false when the
simulation time reaches or exceeds the next output time.

This meant that when the internal timestep exceeded the output
step, `looping` was almost always `false` and so the internal
timestep would saturate.

With this change a timestep adjustment will be considered after
every successful internal timestep.

Note: There may have been a reason why timestep adjustments were
not considered, but I can't remember what it was and it was not
documented.
---
 src/solver/impls/snes/snes.cxx | 93 ++++++++++++++++------------------
 1 file changed, 45 insertions(+), 48 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index de6c54388d..c10095af53 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1069,62 +1069,59 @@ int SNESSolver::run() {
       }
 #endif // PETSC_VERSION_GE(3,20,0)
 
-      if (looping) {
-
-        if (pid_controller) {
-          // Changing the timestep.
-          // Note: The preconditioner depends on the timestep,
-          // so we recalculate the jacobian and the preconditioner
-          //  every time the timestep changes
-
-          timestep = pid(timestep, nl_its);
-
-          // NOTE(malamast): Do we really need this?
-          // Recompute Jacobian (for now)
-          if (saved_jacobian_lag == 0) {
-            SNESGetLagJacobian(snes, &saved_jacobian_lag);
-            SNESSetLagJacobian(snes, 1);
-          }
+      if (pid_controller) {
+        // Changing the timestep.
+        // Note: The preconditioner depends on the timestep,
+        // so we recalculate the jacobian and the preconditioner
+        //  every time the timestep changes
 
-        } else {
+        timestep = pid(timestep, nl_its);
 
-          // Consider changing the timestep.
-          // Note: The preconditioner depends on the timestep,
-          // so if it is not recalculated the it will be less
-          // effective.
-          if ((nl_its <= lower_its) && (timestep < max_timestep)
-              && (steps_since_snes_failure > 2)) {
-            // Increase timestep slightly
-            timestep *= timestep_factor_on_lower_its;
-
-            timestep = std::min(timestep, max_timestep);
-
-            // Note: Setting the SNESJacobianFn to NULL retains
-            // previously set evaluation function.
-            //
-            // The SNES Jacobian is a combination of the RHS Jacobian
-            // and a factor involving the timestep.
-            // Depends on equation_form
-            // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring);
-
-            if (static_cast<BoutReal>(lin_its) / nl_its > 4) {
-              // Recompute Jacobian (for now)
-              if (saved_jacobian_lag == 0) {
-                SNESGetLagJacobian(snes, &saved_jacobian_lag);
-                SNESSetLagJacobian(snes, 1);
-              }
-            }
+        // NOTE(malamast): Do we really need this?
+        // Recompute Jacobian (for now)
+        if (saved_jacobian_lag == 0) {
+          SNESGetLagJacobian(snes, &saved_jacobian_lag);
+          SNESSetLagJacobian(snes, 1);
+        }
 
-          } else if (nl_its >= upper_its) {
-            // Reduce timestep slightly
-            timestep *= timestep_factor_on_upper_its;
+      } else {
 
-            // Recompute Jacobian
+        // Consider changing the timestep.
+        // Note: The preconditioner depends on the timestep,
+        // so if it is not recalculated the it will be less
+        // effective.
+        if ((nl_its <= lower_its) && (timestep < max_timestep)
+            && (steps_since_snes_failure > 2)) {
+          // Increase timestep slightly
+          timestep *= timestep_factor_on_lower_its;
+
+          timestep = std::min(timestep, max_timestep);
+
+          // Note: Setting the SNESJacobianFn to NULL retains
+          // previously set evaluation function.
+          //
+          // The SNES Jacobian is a combination of the RHS Jacobian
+          // and a factor involving the timestep.
+          // Depends on equation_form
+          // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring);
+
+          if (static_cast<BoutReal>(lin_its) / nl_its > 4) {
+            // Recompute Jacobian (for now)
             if (saved_jacobian_lag == 0) {
               SNESGetLagJacobian(snes, &saved_jacobian_lag);
               SNESSetLagJacobian(snes, 1);
             }
           }
+
+        } else if (nl_its >= upper_its) {
+          // Reduce timestep slightly
+          timestep *= timestep_factor_on_upper_its;
+
+          // Recompute Jacobian
+          if (saved_jacobian_lag == 0) {
+            SNESGetLagJacobian(snes, &saved_jacobian_lag);
+            SNESSetLagJacobian(snes, 1);
+          }
         }
       }
       snes_failures = 0;
@@ -1132,7 +1129,7 @@ int SNESSolver::run() {
 
     if (!matrix_free) {
       ASSERT2(simtime >= target);
-      ASSERT2(simtime - dt < target);
+      ASSERT2(simtime - dt <= target);
       // Stepped over output timestep => Interpolate
       // snes_x is the solution at t = simtime
       // x0 is the solution at t = simtime - dt