From 5f0720a5b7831568089118c9bba75ee20006a2c6 Mon Sep 17 00:00:00 2001
From: malamast <malamas.tsagkaridis@gmail.com>
Date: Wed, 2 Jul 2025 18:15:03 -0700
Subject: [PATCH 01/15] SNES solver: Added a PID controller to update the
 timestep of the backward euler SNES solver.

The user sets target_its a desired number of nonlinear iterations and the timestep is updated based on the ratio of that number and the actual number on nonlinear iteration of each Newton step. The exponent can be tuned to the problem needs. A smaller exponent results in smaller changes of dt.
---
 src/solver/impls/snes/snes.cxx | 144 +++++++++++++++++++++++++--------
 src/solver/impls/snes/snes.hxx |  14 ++++
 2 files changed, 125 insertions(+), 33 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 15099f9e2a..8b0b09fea2 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -118,6 +118,9 @@ SNESSolver::SNESSolver(Options* opts)
       maxits((*options)["max_nonlinear_iterations"]
                  .doc("Maximum number of nonlinear iterations per SNES solve")
                  .withDefault(50)),
+      maxf((*options)["maxf"]
+                 .doc("Maximum number of function evaluations per SNES solve")
+                 .withDefault(10000)),                 
       lower_its((*options)["lower_its"]
                     .doc("Iterations below which the next timestep is increased")
                     .withDefault(static_cast<int>(maxits * 0.5))),
@@ -133,6 +136,17 @@ SNESSolver::SNESSolver(Options* opts)
       timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"]
                                    .doc("Multiply timestep if iterations are below lower_its")
                                    .withDefault(1.4)),
+      pidController(
+          (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
+      target_its((*options)["target_its"]
+                    .doc("Target snes iterations")
+                    .withDefault(static_cast<int>(7))),
+      kP((*options)["kP"].doc("Proportional PID parameter")
+                         .withDefault(0.7)),
+      kI((*options)["kI"].doc("Integral PID parameter")
+                         .withDefault(0.3)),
+      kD((*options)["kD"].doc("Derivative PID parameter")
+                         .withDefault(0.2)),   
       diagnose(
           (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)),
       diagnose_failures((*options)["diagnose_failures"]
@@ -632,14 +646,21 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    // Set Jacobian and preconditioner to persist across time steps
-    SNESSetLagJacobianPersists(snes, PETSC_TRUE);
-    SNESSetLagPreconditionerPersists(snes, PETSC_TRUE);
+    if (pidController){
+      nl_its_prev = target_its;
+      nl_its_prev2 = target_its;
+      SNESSetLagJacobianPersists(snes, PETSC_FALSE);
+      SNESSetLagPreconditionerPersists(snes, PETSC_FALSE);
+    } else {
+      // Set Jacobian and preconditioner to persist across time steps
+      SNESSetLagJacobianPersists(snes, PETSC_TRUE);
+      SNESSetLagPreconditionerPersists(snes, PETSC_TRUE);
+    }
     SNESSetLagPreconditioner(snes, 1); // Rebuild when Jacobian is rebuilt
   }
 
   // Set tolerances
-  SNESSetTolerances(snes, atol, rtol, stol, maxits, PETSC_DEFAULT);
+  SNESSetTolerances(snes, atol, rtol, stol, maxits, maxf);
 
   // Force SNES to take at least one nonlinear iteration.
   // This may prevent the solver from getting stuck in false steady state conditions
@@ -843,6 +864,11 @@ int SNESSolver::run() {
         VecAXPBY(snes_x, -beta, (1. + beta), x1);
       }
 
+
+      if (pidController){
+        SNESSetLagJacobian(snes, lag_jacobian);
+      }
+
       // Run the solver
       PetscErrorCode ierr = SNESSolve(snes, nullptr, snes_x);
 
@@ -900,6 +926,7 @@ int SNESSolver::run() {
           updateColoring();
           jacobian_pruned = false; // Reset flag. Will be set after pruning.
         }
+
         if (saved_jacobian_lag == 0) {
           // This triggers a Jacobian recalculation
           SNESGetLagJacobian(snes, &saved_jacobian_lag);
@@ -1033,45 +1060,67 @@ int SNESSolver::run() {
 #endif // PETSC_VERSION_GE(3,20,0)
 
       if (looping) {
-        // 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;
-
-          if (timestep > max_timestep) {
-            timestep = max_timestep;
+
+        if (pidController){
+          // 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);
           }
 
-          // 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)
+        } else {
+          
+          // 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;
+
+            if (timestep > max_timestep) {
+              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);
             }
           }
 
-        } 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;
     } while (looping);
@@ -1386,4 +1435,33 @@ void SNESSolver::updateColoring() {
   }
 }
 
+
+BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
+
+  /* ---------- multiplicative PID factors ---------- */
+  BoutReal facP = std::pow(double(target_its) / double(nl_its),  kP);
+  BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
+  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)  / double(nl_its) / double(nl_its_prev2),  kD);
+
+  // clamp groth factor to avoid huge changes
+  BoutReal fac = facP * facI * facD;
+  if (fac < 0.2)
+    fac = 0.2;
+  else if (fac > 5.0)
+    fac = 5.0;
+
+  /* ---------- update timestep and history ---------- */
+  BoutReal dt_new = timestep * fac;
+
+  if (dt_new > max_timestep) {
+    dt_new = max_timestep;
+  }
+
+  nl_its_prev2 = nl_its_prev;
+  nl_its_prev  = static_cast<int>(nl_its);
+
+  return dt_new;
+}
+
+
 #endif // BOUT_HAS_PETSC
diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx
index cba2352582..02e2368d71 100644
--- a/src/solver/impls/snes/snes.hxx
+++ b/src/solver/impls/snes/snes.hxx
@@ -98,6 +98,7 @@ private:
   BoutReal atol; ///< Absolute tolerance
   BoutReal rtol; ///< Relative tolerance
   BoutReal stol; ///< Convergence tolerance
+  int maxf;                 ///< Maximum number of function evaluations allowed in the solver (default: 10000)
 
   int maxits;               ///< Maximum nonlinear iterations
   int lower_its, upper_its; ///< Limits on iterations for timestep adjustment
@@ -106,6 +107,19 @@ private:
   BoutReal timestep_factor_on_upper_its;
   BoutReal timestep_factor_on_lower_its;
 
+  ///< PID controller parameters
+  bool pidController;  ///< Use PID controller?
+  int target_its;  ///< Target number of nonlinear iterations for the PID controller.
+  ///< Use with caution! Not tested values. 
+  BoutReal kP; ///< (0.6 - 0.8) Proportional parameter (main response to current step)
+  BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes)
+  BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional) 
+
+  int nl_its_prev;
+  int nl_its_prev2;
+
+  BoutReal pid(BoutReal timestep, int nl_its); ///< Updates the timestep
+
   bool diagnose;          ///< Output additional diagnostics
   bool diagnose_failures; ///< Print diagnostics on SNES failures
 

From 656dbdcf6120810c5abc3e24a7e6ba507e398c12 Mon Sep 17 00:00:00 2001
From: malamast <42969991+malamast@users.noreply.github.com>
Date: Thu, 3 Jul 2025 02:28:18 +0000
Subject: [PATCH 02/15] Apply clang-format changes

---
 src/solver/impls/snes/snes.cxx | 57 ++++++++++++++++------------------
 src/solver/impls/snes/snes.hxx | 10 +++---
 2 files changed, 31 insertions(+), 36 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 8b0b09fea2..8b763acf13 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -119,8 +119,8 @@ SNESSolver::SNESSolver(Options* opts)
                  .doc("Maximum number of nonlinear iterations per SNES solve")
                  .withDefault(50)),
       maxf((*options)["maxf"]
-                 .doc("Maximum number of function evaluations per SNES solve")
-                 .withDefault(10000)),                 
+               .doc("Maximum number of function evaluations per SNES solve")
+               .withDefault(10000)),
       lower_its((*options)["lower_its"]
                     .doc("Iterations below which the next timestep is increased")
                     .withDefault(static_cast<int>(maxits * 0.5))),
@@ -128,25 +128,24 @@ SNESSolver::SNESSolver(Options* opts)
                     .doc("Iterations above which the next timestep is reduced")
                     .withDefault(static_cast<int>(maxits * 0.8))),
       timestep_factor_on_failure((*options)["timestep_factor_on_failure"]
-                                 .doc("Multiply timestep on convergence failure")
-                                 .withDefault(0.5)),
-      timestep_factor_on_upper_its((*options)["timestep_factor_on_upper_its"]
-                                   .doc("Multiply timestep if iterations exceed upper_its")
-                                   .withDefault(0.9)),
-      timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"]
-                                   .doc("Multiply timestep if iterations are below lower_its")
-                                   .withDefault(1.4)),
+                                     .doc("Multiply timestep on convergence failure")
+                                     .withDefault(0.5)),
+      timestep_factor_on_upper_its(
+          (*options)["timestep_factor_on_upper_its"]
+              .doc("Multiply timestep if iterations exceed upper_its")
+              .withDefault(0.9)),
+      timestep_factor_on_lower_its(
+          (*options)["timestep_factor_on_lower_its"]
+              .doc("Multiply timestep if iterations are below lower_its")
+              .withDefault(1.4)),
       pidController(
           (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
       target_its((*options)["target_its"]
-                    .doc("Target snes iterations")
-                    .withDefault(static_cast<int>(7))),
-      kP((*options)["kP"].doc("Proportional PID parameter")
-                         .withDefault(0.7)),
-      kI((*options)["kI"].doc("Integral PID parameter")
-                         .withDefault(0.3)),
-      kD((*options)["kD"].doc("Derivative PID parameter")
-                         .withDefault(0.2)),   
+                     .doc("Target snes iterations")
+                     .withDefault(static_cast<int>(7))),
+      kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
+      kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
+      kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),
       diagnose(
           (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)),
       diagnose_failures((*options)["diagnose_failures"]
@@ -646,7 +645,7 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    if (pidController){
+    if (pidController) {
       nl_its_prev = target_its;
       nl_its_prev2 = target_its;
       SNESSetLagJacobianPersists(snes, PETSC_FALSE);
@@ -864,8 +863,7 @@ int SNESSolver::run() {
         VecAXPBY(snes_x, -beta, (1. + beta), x1);
       }
 
-
-      if (pidController){
+      if (pidController) {
         SNESSetLagJacobian(snes, lag_jacobian);
       }
 
@@ -1061,7 +1059,7 @@ int SNESSolver::run() {
 
       if (looping) {
 
-        if (pidController){
+        if (pidController) {
           // Changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so we recalculate the jacobian and the preconditioner
@@ -1076,9 +1074,8 @@ int SNESSolver::run() {
             SNESSetLagJacobian(snes, 1);
           }
 
-
         } else {
-          
+
           // Consider changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so if it is not recalculated the it will be less
@@ -1118,9 +1115,7 @@ int SNESSolver::run() {
               SNESSetLagJacobian(snes, 1);
             }
           }
-
         }
-
       }
       snes_failures = 0;
     } while (looping);
@@ -1435,13 +1430,14 @@ void SNESSolver::updateColoring() {
   }
 }
 
-
 BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
 
   /* ---------- multiplicative PID factors ---------- */
-  BoutReal facP = std::pow(double(target_its) / double(nl_its),  kP);
+  BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
   BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)  / double(nl_its) / double(nl_its_prev2),  kD);
+  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
+                               / double(nl_its_prev2),
+                           kD);
 
   // clamp groth factor to avoid huge changes
   BoutReal fac = facP * facI * facD;
@@ -1458,10 +1454,9 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   }
 
   nl_its_prev2 = nl_its_prev;
-  nl_its_prev  = static_cast<int>(nl_its);
+  nl_its_prev = static_cast<int>(nl_its);
 
   return dt_new;
 }
 
-
 #endif // BOUT_HAS_PETSC
diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx
index 02e2368d71..8a0b5bcecb 100644
--- a/src/solver/impls/snes/snes.hxx
+++ b/src/solver/impls/snes/snes.hxx
@@ -98,7 +98,7 @@ private:
   BoutReal atol; ///< Absolute tolerance
   BoutReal rtol; ///< Relative tolerance
   BoutReal stol; ///< Convergence tolerance
-  int maxf;                 ///< Maximum number of function evaluations allowed in the solver (default: 10000)
+  int maxf; ///< Maximum number of function evaluations allowed in the solver (default: 10000)
 
   int maxits;               ///< Maximum nonlinear iterations
   int lower_its, upper_its; ///< Limits on iterations for timestep adjustment
@@ -108,12 +108,12 @@ private:
   BoutReal timestep_factor_on_lower_its;
 
   ///< PID controller parameters
-  bool pidController;  ///< Use PID controller?
-  int target_its;  ///< Target number of nonlinear iterations for the PID controller.
-  ///< Use with caution! Not tested values. 
+  bool pidController; ///< Use PID controller?
+  int target_its;     ///< Target number of nonlinear iterations for the PID controller.
+  ///< Use with caution! Not tested values.
   BoutReal kP; ///< (0.6 - 0.8) Proportional parameter (main response to current step)
   BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes)
-  BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional) 
+  BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional)
 
   int nl_its_prev;
   int nl_its_prev2;

From 5a240c85eb6fa8bc22b0351342a50e7205517d4d Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Mon, 7 Jul 2025 13:10:51 -0700
Subject: [PATCH 03/15] Update src/solver/impls/snes/snes.cxx

Code changes suggested by ZedThree

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 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 8b763acf13..678713f445 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -142,7 +142,7 @@ SNESSolver::SNESSolver(Options* opts)
           (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
       target_its((*options)["target_its"]
                      .doc("Target snes iterations")
-                     .withDefault(static_cast<int>(7))),
+                     .withDefault(7)),
       kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
       kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
       kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),

From 9276da0e3fa50b7a9f613dabf011b018fa16004e Mon Sep 17 00:00:00 2001
From: malamast <42969991+malamast@users.noreply.github.com>
Date: Mon, 7 Jul 2025 20:12:06 +0000
Subject: [PATCH 04/15] Apply clang-format changes

---
 src/bout++.cxx                 | 2 +-
 src/solver/impls/snes/snes.cxx | 4 +---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/bout++.cxx b/src/bout++.cxx
index 7f23cf5f91..4e9e37fe8d 100644
--- a/src/bout++.cxx
+++ b/src/bout++.cxx
@@ -84,7 +84,7 @@ const char DEFAULT_DIR[] = "data";
 // Define S_ISDIR if not defined by system headers (that is, MSVC)
 // Taken from https://github.com/curl/curl/blob/e59540139a398dc70fde6aec487b19c5085105af/lib/curl_setup.h#L748-L751
 #if !defined(S_ISDIR) && defined(S_IFMT) && defined(S_IFDIR)
-#define S_ISDIR(m) (((m)&S_IFMT) == S_IFDIR)
+#define S_ISDIR(m) (((m) & S_IFMT) == S_IFDIR)
 #endif
 
 #ifdef _MSC_VER
diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 678713f445..9906274a6e 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -140,9 +140,7 @@ SNESSolver::SNESSolver(Options* opts)
               .withDefault(1.4)),
       pidController(
           (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
-      target_its((*options)["target_its"]
-                     .doc("Target snes iterations")
-                     .withDefault(7)),
+      target_its((*options)["target_its"].doc("Target snes iterations").withDefault(7)),
       kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
       kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
       kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),

From 210a6bc51499aaf39a01ee46ce4408b5909692da Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 11:34:07 -0700
Subject: [PATCH 05/15] Update src/solver/impls/snes/snes.cxx

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 src/solver/impls/snes/snes.cxx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 9906274a6e..387c908cb8 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1431,9 +1431,9 @@ void SNESSolver::updateColoring() {
 BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
 
   /* ---------- multiplicative PID factors ---------- */
-  BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
-  BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
+  const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
+  const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
+  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
                                / double(nl_its_prev2),
                            kD);
 

From e881832ad7d0d5447886b0787af8ab179d9c833d Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 11:34:26 -0700
Subject: [PATCH 06/15] Update src/solver/impls/snes/snes.cxx

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 src/solver/impls/snes/snes.cxx | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 387c908cb8..a0f9c06801 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1437,12 +1437,8 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
                                / double(nl_its_prev2),
                            kD);
 
-  // clamp groth factor to avoid huge changes
-  BoutReal fac = facP * facI * facD;
-  if (fac < 0.2)
-    fac = 0.2;
-  else if (fac > 5.0)
-    fac = 5.0;
+  // clamp growth factor to avoid huge changes
+  const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);
 
   /* ---------- update timestep and history ---------- */
   BoutReal dt_new = timestep * fac;

From 133980416d2b266e15a718b3ef045b7b355a37e4 Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 11:34:37 -0700
Subject: [PATCH 07/15] Update src/solver/impls/snes/snes.cxx

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 src/solver/impls/snes/snes.cxx | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index a0f9c06801..9d9a396c3c 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1441,11 +1441,7 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);
 
   /* ---------- update timestep and history ---------- */
-  BoutReal dt_new = timestep * fac;
-
-  if (dt_new > max_timestep) {
-    dt_new = max_timestep;
-  }
+  const BoutReal dt_new = std::min(timestep * fac, max_timestep);
 
   nl_its_prev2 = nl_its_prev;
   nl_its_prev = static_cast<int>(nl_its);

From 687d21c5da8e40c68ff5d7306ae769a058ce9038 Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 11:34:47 -0700
Subject: [PATCH 08/15] Update src/solver/impls/snes/snes.cxx

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 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 9d9a396c3c..89910f1655 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1444,7 +1444,7 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   const BoutReal dt_new = std::min(timestep * fac, max_timestep);
 
   nl_its_prev2 = nl_its_prev;
-  nl_its_prev = static_cast<int>(nl_its);
+  nl_its_prev = nl_its;
 
   return dt_new;
 }

From 0c8ef879a944b2d58fdbdc04ac74f69f2c37c33f Mon Sep 17 00:00:00 2001
From: Malamas Tsagkaridis <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 11:34:57 -0700
Subject: [PATCH 09/15] Update src/solver/impls/snes/snes.cxx

Co-authored-by: Peter Hill <peter.hill@york.ac.uk>
---
 src/solver/impls/snes/snes.cxx | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 89910f1655..0b61320f34 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1083,9 +1083,7 @@ int SNESSolver::run() {
             // Increase timestep slightly
             timestep *= timestep_factor_on_lower_its;
 
-            if (timestep > max_timestep) {
-              timestep = max_timestep;
-            }
+            timestep = std::min(timestep, max_timestep);
 
             // Note: Setting the SNESJacobianFn to NULL retains
             // previously set evaluation function.

From fd6f18943cf2d715a507ce50cb7ff702a21b6196 Mon Sep 17 00:00:00 2001
From: malamast <42969991+malamast@users.noreply.github.com>
Date: Tue, 8 Jul 2025 18:36:20 +0000
Subject: [PATCH 10/15] Apply clang-format changes

---
 src/solver/impls/snes/snes.cxx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 0b61320f34..e14d3d4569 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1431,9 +1431,9 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   /* ---------- multiplicative PID factors ---------- */
   const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
   const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
-                               / double(nl_its_prev2),
-                           kD);
+  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)
+                                     / double(nl_its) / double(nl_its_prev2),
+                                 kD);
 
   // clamp growth factor to avoid huge changes
   const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);

From 8619876339753cc0c5f8c1b4810adb274a57fc26 Mon Sep 17 00:00:00 2001
From: malamast <malamas.tsagkaridis@gmail.com>
Date: Thu, 17 Jul 2025 09:56:23 -0700
Subject: [PATCH 11/15] snes: I addied some changes made on my remote branch by
 clang-tidy. I added the line
 MatSetOption(Jfd,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE); after the creation of
 Jfd to tell petsc to keep the pattern of the non zero elements in memory.
 Also, I changed the default coloring type to greedy which is more suited for
 large parallel problems.

---
 src/solver/impls/snes/snes.cxx | 82 +++++++++++++++++++++-------------
 1 file changed, 51 insertions(+), 31 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index e14d3d4569..4f6c11fb42 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -119,8 +119,8 @@ SNESSolver::SNESSolver(Options* opts)
                  .doc("Maximum number of nonlinear iterations per SNES solve")
                  .withDefault(50)),
       maxf((*options)["maxf"]
-               .doc("Maximum number of function evaluations per SNES solve")
-               .withDefault(10000)),
+                 .doc("Maximum number of function evaluations per SNES solve")
+                 .withDefault(10000)),                 
       lower_its((*options)["lower_its"]
                     .doc("Iterations below which the next timestep is increased")
                     .withDefault(static_cast<int>(maxits * 0.5))),
@@ -128,22 +128,25 @@ SNESSolver::SNESSolver(Options* opts)
                     .doc("Iterations above which the next timestep is reduced")
                     .withDefault(static_cast<int>(maxits * 0.8))),
       timestep_factor_on_failure((*options)["timestep_factor_on_failure"]
-                                     .doc("Multiply timestep on convergence failure")
-                                     .withDefault(0.5)),
-      timestep_factor_on_upper_its(
-          (*options)["timestep_factor_on_upper_its"]
-              .doc("Multiply timestep if iterations exceed upper_its")
-              .withDefault(0.9)),
-      timestep_factor_on_lower_its(
-          (*options)["timestep_factor_on_lower_its"]
-              .doc("Multiply timestep if iterations are below lower_its")
-              .withDefault(1.4)),
+                                 .doc("Multiply timestep on convergence failure")
+                                 .withDefault(0.5)),
+      timestep_factor_on_upper_its((*options)["timestep_factor_on_upper_its"]
+                                   .doc("Multiply timestep if iterations exceed upper_its")
+                                   .withDefault(0.9)),
+      timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"]
+                                   .doc("Multiply timestep if iterations are below lower_its")
+                                   .withDefault(1.4)),
       pidController(
           (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
-      target_its((*options)["target_its"].doc("Target snes iterations").withDefault(7)),
-      kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
-      kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
-      kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),
+      target_its((*options)["target_its"]
+                    .doc("Target snes iterations")
+                    .withDefault(static_cast<int>(7))),
+      kP((*options)["kP"].doc("Proportional PID parameter")
+                         .withDefault(0.7)),
+      kI((*options)["kI"].doc("Integral PID parameter")
+                         .withDefault(0.3)),
+      kD((*options)["kD"].doc("Derivative PID parameter")
+                         .withDefault(0.2)),   
       diagnose(
           (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)),
       diagnose_failures((*options)["diagnose_failures"]
@@ -351,6 +354,7 @@ int SNESSolver::init() {
 
       // Set size of Matrix on each processor to nlocal x nlocal
       MatCreate(BoutComm::get(), &Jfd);
+      MatSetOption(Jfd,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
       MatSetSizes(Jfd, nlocal, nlocal, PETSC_DETERMINE, PETSC_DETERMINE);
       MatSetFromOptions(Jfd);
       // Determine which row/columns of the matrix are locally owned
@@ -643,7 +647,7 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    if (pidController) {
+    if (pidController){
       nl_its_prev = target_its;
       nl_its_prev2 = target_its;
       SNESSetLagJacobianPersists(snes, PETSC_FALSE);
@@ -861,7 +865,8 @@ int SNESSolver::run() {
         VecAXPBY(snes_x, -beta, (1. + beta), x1);
       }
 
-      if (pidController) {
+
+      if (pidController){
         SNESSetLagJacobian(snes, lag_jacobian);
       }
 
@@ -1057,7 +1062,7 @@ int SNESSolver::run() {
 
       if (looping) {
 
-        if (pidController) {
+        if (pidController){
           // Changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so we recalculate the jacobian and the preconditioner
@@ -1072,8 +1077,9 @@ int SNESSolver::run() {
             SNESSetLagJacobian(snes, 1);
           }
 
-        } else {
 
+        } else {
+          
           // Consider changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so if it is not recalculated the it will be less
@@ -1083,7 +1089,9 @@ int SNESSolver::run() {
             // Increase timestep slightly
             timestep *= timestep_factor_on_lower_its;
 
-            timestep = std::min(timestep, max_timestep);
+            if (timestep > max_timestep) {
+              timestep = max_timestep;
+            }
 
             // Note: Setting the SNESJacobianFn to NULL retains
             // previously set evaluation function.
@@ -1111,7 +1119,9 @@ int SNESSolver::run() {
               SNESSetLagJacobian(snes, 1);
             }
           }
+
         }
+
       }
       snes_failures = 0;
     } while (looping);
@@ -1400,7 +1410,9 @@ void SNESSolver::updateColoring() {
   // Re-calculate the coloring
   MatColoring coloring = NULL;
   MatColoringCreate(Jfd, &coloring);
-  MatColoringSetType(coloring, MATCOLORINGSL);
+  // MatColoringSetType(coloring, MATCOLORINGSL);  // Serial algorithm. Better for smale-to-medium size problems.
+  MatColoringSetType(coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
+  // MatColoringSetType(coloring, MATCOLORINGJP);  // This didn't work
   MatColoringSetFromOptions(coloring);
 
   // Calculate new index sets
@@ -1426,25 +1438,33 @@ void SNESSolver::updateColoring() {
   }
 }
 
+
 BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
 
   /* ---------- multiplicative PID factors ---------- */
-  const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
-  const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)
-                                     / double(nl_its) / double(nl_its_prev2),
-                                 kD);
+  BoutReal facP = std::pow(double(target_its) / double(nl_its),  kP);
+  BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
+  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)  / double(nl_its) / double(nl_its_prev2),  kD);
 
-  // clamp growth factor to avoid huge changes
-  const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);
+  // clamp groth factor to avoid huge changes
+  BoutReal fac = facP * facI * facD;
+  if (fac < 0.2)
+    fac = 0.2;
+  else if (fac > 5.0)
+    fac = 5.0;
 
   /* ---------- update timestep and history ---------- */
-  const BoutReal dt_new = std::min(timestep * fac, max_timestep);
+  BoutReal dt_new = timestep * fac;
+
+  if (dt_new > max_timestep) {
+    dt_new = max_timestep;
+  }
 
   nl_its_prev2 = nl_its_prev;
-  nl_its_prev = nl_its;
+  nl_its_prev  = static_cast<int>(nl_its);
 
   return dt_new;
 }
 
+
 #endif // BOUT_HAS_PETSC

From c20863993e7e87565cecef5a16eb18a97ab9518a Mon Sep 17 00:00:00 2001
From: malamast <42969991+malamast@users.noreply.github.com>
Date: Thu, 17 Jul 2025 17:59:07 +0000
Subject: [PATCH 12/15] Apply clang-format changes

---
 src/solver/impls/snes/snes.cxx | 62 ++++++++++++++++------------------
 1 file changed, 29 insertions(+), 33 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 4f6c11fb42..0b57a52d39 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -119,8 +119,8 @@ SNESSolver::SNESSolver(Options* opts)
                  .doc("Maximum number of nonlinear iterations per SNES solve")
                  .withDefault(50)),
       maxf((*options)["maxf"]
-                 .doc("Maximum number of function evaluations per SNES solve")
-                 .withDefault(10000)),                 
+               .doc("Maximum number of function evaluations per SNES solve")
+               .withDefault(10000)),
       lower_its((*options)["lower_its"]
                     .doc("Iterations below which the next timestep is increased")
                     .withDefault(static_cast<int>(maxits * 0.5))),
@@ -128,25 +128,24 @@ SNESSolver::SNESSolver(Options* opts)
                     .doc("Iterations above which the next timestep is reduced")
                     .withDefault(static_cast<int>(maxits * 0.8))),
       timestep_factor_on_failure((*options)["timestep_factor_on_failure"]
-                                 .doc("Multiply timestep on convergence failure")
-                                 .withDefault(0.5)),
-      timestep_factor_on_upper_its((*options)["timestep_factor_on_upper_its"]
-                                   .doc("Multiply timestep if iterations exceed upper_its")
-                                   .withDefault(0.9)),
-      timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"]
-                                   .doc("Multiply timestep if iterations are below lower_its")
-                                   .withDefault(1.4)),
+                                     .doc("Multiply timestep on convergence failure")
+                                     .withDefault(0.5)),
+      timestep_factor_on_upper_its(
+          (*options)["timestep_factor_on_upper_its"]
+              .doc("Multiply timestep if iterations exceed upper_its")
+              .withDefault(0.9)),
+      timestep_factor_on_lower_its(
+          (*options)["timestep_factor_on_lower_its"]
+              .doc("Multiply timestep if iterations are below lower_its")
+              .withDefault(1.4)),
       pidController(
           (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
       target_its((*options)["target_its"]
-                    .doc("Target snes iterations")
-                    .withDefault(static_cast<int>(7))),
-      kP((*options)["kP"].doc("Proportional PID parameter")
-                         .withDefault(0.7)),
-      kI((*options)["kI"].doc("Integral PID parameter")
-                         .withDefault(0.3)),
-      kD((*options)["kD"].doc("Derivative PID parameter")
-                         .withDefault(0.2)),   
+                     .doc("Target snes iterations")
+                     .withDefault(static_cast<int>(7))),
+      kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
+      kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
+      kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),
       diagnose(
           (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)),
       diagnose_failures((*options)["diagnose_failures"]
@@ -354,7 +353,7 @@ int SNESSolver::init() {
 
       // Set size of Matrix on each processor to nlocal x nlocal
       MatCreate(BoutComm::get(), &Jfd);
-      MatSetOption(Jfd,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
+      MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
       MatSetSizes(Jfd, nlocal, nlocal, PETSC_DETERMINE, PETSC_DETERMINE);
       MatSetFromOptions(Jfd);
       // Determine which row/columns of the matrix are locally owned
@@ -647,7 +646,7 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    if (pidController){
+    if (pidController) {
       nl_its_prev = target_its;
       nl_its_prev2 = target_its;
       SNESSetLagJacobianPersists(snes, PETSC_FALSE);
@@ -865,8 +864,7 @@ int SNESSolver::run() {
         VecAXPBY(snes_x, -beta, (1. + beta), x1);
       }
 
-
-      if (pidController){
+      if (pidController) {
         SNESSetLagJacobian(snes, lag_jacobian);
       }
 
@@ -1062,7 +1060,7 @@ int SNESSolver::run() {
 
       if (looping) {
 
-        if (pidController){
+        if (pidController) {
           // Changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so we recalculate the jacobian and the preconditioner
@@ -1077,9 +1075,8 @@ int SNESSolver::run() {
             SNESSetLagJacobian(snes, 1);
           }
 
-
         } else {
-          
+
           // Consider changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so if it is not recalculated the it will be less
@@ -1119,9 +1116,7 @@ int SNESSolver::run() {
               SNESSetLagJacobian(snes, 1);
             }
           }
-
         }
-
       }
       snes_failures = 0;
     } while (looping);
@@ -1411,7 +1406,8 @@ void SNESSolver::updateColoring() {
   MatColoring coloring = NULL;
   MatColoringCreate(Jfd, &coloring);
   // MatColoringSetType(coloring, MATCOLORINGSL);  // Serial algorithm. Better for smale-to-medium size problems.
-  MatColoringSetType(coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
+  MatColoringSetType(
+      coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
   // MatColoringSetType(coloring, MATCOLORINGJP);  // This didn't work
   MatColoringSetFromOptions(coloring);
 
@@ -1438,13 +1434,14 @@ void SNESSolver::updateColoring() {
   }
 }
 
-
 BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
 
   /* ---------- multiplicative PID factors ---------- */
-  BoutReal facP = std::pow(double(target_its) / double(nl_its),  kP);
+  BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
   BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)  / double(nl_its) / double(nl_its_prev2),  kD);
+  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
+                               / double(nl_its_prev2),
+                           kD);
 
   // clamp groth factor to avoid huge changes
   BoutReal fac = facP * facI * facD;
@@ -1461,10 +1458,9 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   }
 
   nl_its_prev2 = nl_its_prev;
-  nl_its_prev  = static_cast<int>(nl_its);
+  nl_its_prev = static_cast<int>(nl_its);
 
   return dt_new;
 }
 
-
 #endif // BOUT_HAS_PETSC

From d816d7819e128925fe806052aeebb5ecdac1eaf2 Mon Sep 17 00:00:00 2001
From: malamast <malamas.tsagkaridis@gmail.com>
Date: Fri, 18 Jul 2025 10:23:33 -0700
Subject: [PATCH 13/15] SNES: changes the pidController flag to pid_controller
 to match the variable naming convention.

---
 src/solver/impls/snes/snes.cxx | 10 +++++-----
 src/solver/impls/snes/snes.hxx |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 4f6c11fb42..8c10052e33 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -136,8 +136,8 @@ SNESSolver::SNESSolver(Options* opts)
       timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"]
                                    .doc("Multiply timestep if iterations are below lower_its")
                                    .withDefault(1.4)),
-      pidController(
-          (*options)["pidController"].doc("Use PID controller?").withDefault(false)),
+      pid_controller(
+          (*options)["pid_controller"].doc("Use PID controller?").withDefault(false)),
       target_its((*options)["target_its"]
                     .doc("Target snes iterations")
                     .withDefault(static_cast<int>(7))),
@@ -647,7 +647,7 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    if (pidController){
+    if (pid_controller){
       nl_its_prev = target_its;
       nl_its_prev2 = target_its;
       SNESSetLagJacobianPersists(snes, PETSC_FALSE);
@@ -866,7 +866,7 @@ int SNESSolver::run() {
       }
 
 
-      if (pidController){
+      if (pid_controller){
         SNESSetLagJacobian(snes, lag_jacobian);
       }
 
@@ -1062,7 +1062,7 @@ int SNESSolver::run() {
 
       if (looping) {
 
-        if (pidController){
+        if (pid_controller){
           // Changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so we recalculate the jacobian and the preconditioner
diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx
index 8a0b5bcecb..e31ca44625 100644
--- a/src/solver/impls/snes/snes.hxx
+++ b/src/solver/impls/snes/snes.hxx
@@ -108,7 +108,7 @@ private:
   BoutReal timestep_factor_on_lower_its;
 
   ///< PID controller parameters
-  bool pidController; ///< Use PID controller?
+  bool pid_controller; ///< Use PID controller?
   int target_its;     ///< Target number of nonlinear iterations for the PID controller.
   ///< Use with caution! Not tested values.
   BoutReal kP; ///< (0.6 - 0.8) Proportional parameter (main response to current step)

From 751a22f6ce2a75047c3e61047bfd4aec5ce62ee2 Mon Sep 17 00:00:00 2001
From: malamast <malamas.tsagkaridis@gmail.com>
Date: Thu, 24 Jul 2025 11:34:02 -0700
Subject: [PATCH 14/15] SNES: reapplied some changes from previous commits that
 accidentally got deleted.

---
 src/solver/impls/snes/snes.cxx | 37 +++++++++++-----------------------
 src/solver/impls/snes/snes.hxx |  2 +-
 2 files changed, 13 insertions(+), 26 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 6b0438bb8b..eda4e9e783 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -140,9 +140,7 @@ SNESSolver::SNESSolver(Options* opts)
               .withDefault(1.4)),
       pid_controller(
           (*options)["pid_controller"].doc("Use PID controller?").withDefault(false)),
-      target_its((*options)["target_its"]
-                     .doc("Target snes iterations")
-                     .withDefault(static_cast<int>(7))),
+      target_its((*options)["target_its"].doc("Target snes iterations").withDefault(7)),
       kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)),
       kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)),
       kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)),
@@ -646,7 +644,7 @@ int SNESSolver::init() {
     // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian
     //       always updates its reference 'u' vector every nonlinear iteration
     SNESSetLagJacobian(snes, lag_jacobian);
-    if (pid_controller){
+    if (pid_controller) {
       nl_its_prev = target_its;
       nl_its_prev2 = target_its;
       SNESSetLagJacobianPersists(snes, PETSC_FALSE);
@@ -865,7 +863,7 @@ int SNESSolver::run() {
       }
 
 
-      if (pid_controller){
+      if (pid_controller) {
         SNESSetLagJacobian(snes, lag_jacobian);
       }
 
@@ -1061,7 +1059,7 @@ int SNESSolver::run() {
 
       if (looping) {
 
-        if (pid_controller){
+        if (pid_controller) {
           // Changing the timestep.
           // Note: The preconditioner depends on the timestep,
           // so we recalculate the jacobian and the preconditioner
@@ -1087,9 +1085,7 @@ int SNESSolver::run() {
             // Increase timestep slightly
             timestep *= timestep_factor_on_lower_its;
 
-            if (timestep > max_timestep) {
-              timestep = max_timestep;
-            }
+            timestep = std::min(timestep, max_timestep);
 
             // Note: Setting the SNESJacobianFn to NULL retains
             // previously set evaluation function.
@@ -1407,8 +1403,7 @@ void SNESSolver::updateColoring() {
   MatColoring coloring = NULL;
   MatColoringCreate(Jfd, &coloring);
   // MatColoringSetType(coloring, MATCOLORINGSL);  // Serial algorithm. Better for smale-to-medium size problems.
-  MatColoringSetType(
-      coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
+  MatColoringSetType(coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
   // MatColoringSetType(coloring, MATCOLORINGJP);  // This didn't work
   MatColoringSetFromOptions(coloring);
 
@@ -1438,28 +1433,20 @@ void SNESSolver::updateColoring() {
 BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
 
   /* ---------- multiplicative PID factors ---------- */
-  BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
-  BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
+  const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
+  const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
+  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
                                / double(nl_its_prev2),
                            kD);
 
   // clamp groth factor to avoid huge changes
-  BoutReal fac = facP * facI * facD;
-  if (fac < 0.2)
-    fac = 0.2;
-  else if (fac > 5.0)
-    fac = 5.0;
+  const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);
 
   /* ---------- update timestep and history ---------- */
-  BoutReal dt_new = timestep * fac;
-
-  if (dt_new > max_timestep) {
-    dt_new = max_timestep;
-  }
+  const BoutReal dt_new = std::min(timestep * fac, max_timestep);
 
   nl_its_prev2 = nl_its_prev;
-  nl_its_prev = static_cast<int>(nl_its);
+  nl_its_prev = nl_its;
 
   return dt_new;
 }
diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx
index e31ca44625..bd942f09ff 100644
--- a/src/solver/impls/snes/snes.hxx
+++ b/src/solver/impls/snes/snes.hxx
@@ -109,7 +109,7 @@ private:
 
   ///< PID controller parameters
   bool pid_controller; ///< Use PID controller?
-  int target_its;     ///< Target number of nonlinear iterations for the PID controller.
+  int target_its;      ///< Target number of nonlinear iterations for the PID controller.
   ///< Use with caution! Not tested values.
   BoutReal kP; ///< (0.6 - 0.8) Proportional parameter (main response to current step)
   BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes)

From 865c0ab26987f393b9e8ab9294603259a501e3f1 Mon Sep 17 00:00:00 2001
From: malamast <42969991+malamast@users.noreply.github.com>
Date: Thu, 24 Jul 2025 19:35:33 +0000
Subject: [PATCH 15/15] Apply clang-format changes

---
 src/solver/impls/snes/snes.cxx | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index eda4e9e783..72891941b6 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -862,7 +862,6 @@ int SNESSolver::run() {
         VecAXPBY(snes_x, -beta, (1. + beta), x1);
       }
 
-
       if (pid_controller) {
         SNESSetLagJacobian(snes, lag_jacobian);
       }
@@ -1403,7 +1402,8 @@ void SNESSolver::updateColoring() {
   MatColoring coloring = NULL;
   MatColoringCreate(Jfd, &coloring);
   // MatColoringSetType(coloring, MATCOLORINGSL);  // Serial algorithm. Better for smale-to-medium size problems.
-  MatColoringSetType(coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
+  MatColoringSetType(
+      coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs
   // MatColoringSetType(coloring, MATCOLORINGJP);  // This didn't work
   MatColoringSetFromOptions(coloring);
 
@@ -1435,9 +1435,9 @@ BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) {
   /* ---------- multiplicative PID factors ---------- */
   const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
   const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
-  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
-                               / double(nl_its_prev2),
-                           kD);
+  const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev)
+                                     / double(nl_its) / double(nl_its_prev2),
+                                 kD);
 
   // clamp groth factor to avoid huge changes
   const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);