snes solver: Pseudo Transient Continuation method #3181

bendudson · 2025-10-12T18:38:56Z

~~Setting pseudo_time = true now enables Pseudo-timestepping, in which each cell has a separate timestep.~~

Setting equation_form = pseudo_transient now enables Pseudo-timestepping, in which each cell has a separate timestep.

The timestep in each cell is set inversely proportional to the residual (RMS time-derivative of all quantities in cell).

Recommend enabling pid_controller that multiplies all timesteps by the same factor, to control the number of nonlinear iterations.

Output time uses the minimum timestep in the domain, so simulations shouldn't need to run for as long as when all cells are evolved at the same speed.

Includes PR #3180

Settings that seem to work well for the Hermes-3 examples/tokamak-1D/1D-threshold case:

[solver]
type = snes
equation_form = pseudo_transient    # Enable pseudo-time stepping

pid_controller = true   # Use PID controller to scale all timesteps based on solver convergence
target_its = 10              # Target number of nonlinear iterations
kP = 0.233
kI = 0.133
kD = 0.0
pid_consider_failures = true  # PID controller becomes cautious after solver failures

snes_type = newtonls
ksp_type = fgmres
pc_type = lu
max_nonlinear_iterations = 20   # Limit significantly above target_its to limit failures
atol = 1e-7
rtol = 1e-6
stol = 1e-12
maxf = 200
maxl = 20

matrix_free_operator = true   # Use the actual nonlinear function for the matrix-vector product

lag_jacobian = 5   # With pid_controller=true the Jacobian will be recalculated every timestep
timestep = 1.0       # Starting internal timestep

[petsc]
snes_ksp_ew = false # Eisenstat-Walker method for setting linear tolerances
ksp_type = fgmres
pc_type = lu
pc_factor_mat_solver_type = strumpack
snes_linesearch_type = basic

This documentation of ADflow contains some useful discussion: https://mdolab-adflow.readthedocs-hosted.com/en/latest/solvers.html
The methods used in the ADflow solver seems to be very similar to this: Pseudo-Transient Continuation, matrix-free operator, Inexact Newton-Krylov, Eisenstat-Walker.

When using Eisenstat-Walker the number of nonlinear iterations should be increased: Without it I found target_its=4 was ok, but with EW enabled the solver stalled with small timesteps unless I increased target_its to e.g 10.

Testing Hermes-3 examples/tokamak-1D/1D-threshold on 1 core:

Previous settings (Hermes-3 master branch) take 1m 18s to run and have not yet reached steady state
Pseudo-time settings (above) run in 16s and have reached steady state.

For 2D problems try:

[solver]
diagnose = true
type = snes
equation_form = pseudo_transient    # Enable pseudo-time stepping

pid_controller = true   # Use PID controller to scale all timesteps based on solver convergence
target_its = 4         # Target number of nonlinear iterations
kP = 0.233
kI = 0.133
kD = 0.0

snes_type = newtonls
ksp_type = fgmres
pc_type = lu
max_nonlinear_iterations = 20   # Limit significantly above target_its to limit failures
atol = 1e-7
rtol = 1e-6
stol = 1e-12
maxf = 200
maxl = 20

matrix_free_operator = true   # Use the actual nonlinear function for the matrix-vector product

lag_jacobian = 5   # With pid_controller=true the Jacobian will be recalculated every timestep
timestep = 1.0       # Starting internal timestep

[petsc]
snes_ksp_ew = false # Eisenstat-Walker method for setting linear tolerances
ksp_type = fgmres
pc_type = lu
pc_factor_mat_solver_type = strumpack
snes_linesearch_type = basic

Setting `pseudo_time = true` now enables Pseudo-timestepping, in which each cell has a separate timestep. The timestep in each cell is set inversely proportional to the residual (RMS time-derivative of all quantities in cell). Recommend enabling `pid_controller` that multiplies all timesteps by the same factor, to control the number of nonlinear iterations. Output time uses the minimum timestep in the domain, so simulations shouldn't need to run for as long as when all cells are evolved at the same speed.

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 30. Check the log or trigger a new build to see more.

github-actions · 2025-10-12T18:50:55Z

src/solver/impls/snes/snes.cxx

 #include <bout/globals.hxx>
 #include <bout/msg_stack.hxx>
+#include <bout/output_bout_types.hxx>
 #include <bout/petsc_interface.hxx>


warning: included header output_bout_types.hxx is not used directly [misc-include-cleaner]

Suggested change

#include <bout/petsc_interface.hxx>

#include <bout/petsc_interface.hxx>

src/solver/impls/snes/snes.cxx

Rather than having a separate switch, just use equation_form to select the Pseudo-Transient Continuation method.

The rhs_function returns the time derivatives (residuals) that are used in the PTC method. The snes_function calls rhs_function, then modifies the function depending on equation_form.

github-actions

clang-tidy made some suggestions

src/solver/impls/snes/snes.cxx

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

+
+  // Timestep inversely proportional to the residual, clipped to avoid
+  // rapid changes in timestep.
+  return std::min(


warning: no header providing "std::min" is directly included [misc-include-cleaner]

. ^

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

+  // Timestep inversely proportional to the residual, clipped to avoid
+  // rapid changes in timestep.
+  return std::min(
+      {std::max({pseudo_alpha / current_residual, previous_timestep / 1.5, dt_min_reset}),


warning: no header providing "std::max" is directly included [misc-include-cleaner]

( ^

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

-    CHKERRQ(ierr);
+    PetscCall(VecGetArrayRead(scaled_x, &xdata));
+    // const_cast needed due to load_vars API. Not writing to xdata.
+    load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast to remove const qualifier [cppcoreguidelines-pro-type-const-cast]

. ^

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

-    CHKERRQ(ierr);
+    PetscCall(VecGetArrayRead(x, &xdata));
+    // const_cast needed due to load_vars API. Not writing to xdata.
+    load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast to remove const qualifier [cppcoreguidelines-pro-type-const-cast]

. ^

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

+  PetscCall(VecGetArray(f, &fdata));
+  save_derivs(fdata);
+  PetscCall(VecRestoreArray(f, &fdata));
+  return PETSC_SUCCESS;


warning: no header providing "PETSC_SUCCESS" is directly included [misc-include-cleaner]

; ^

github-actions · 2025-10-13T18:24:23Z

src/solver/impls/snes/snes.cxx

  // Scale the Jacobian rows by multiplying on the left by 1/norm
-  ierr = MatDiagonalScale(Jac_new, jac_row_inv_norms, nullptr);
-  CHKERRQ(ierr);
+  PetscCall(MatDiagonalScale(Jac_new, jac_row_inv_norms, nullptr));


warning: no header providing "MatDiagonalScale" is directly included [misc-include-cleaner]

1/norm ^

Switches between different strategies for choosing cell time steps.

Minor tidying

Trying to simplify the `run()` function by splitting out pieces into separate functions.

github-actions

clang-tidy made some suggestions

github-actions · 2025-10-13T22:17:51Z

src/solver/impls/snes/snes.cxx

  SNESSolver* s;
-  int ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: do not use reinterpret_cast [cppcoreguidelines-pro-type-reinterpret-cast]

PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); ^

github-actions · 2025-10-13T22:17:52Z

src/solver/impls/snes/snes.cxx

  SNESSolver* s;
-  int ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: multilevel pointer conversion from 'void **' to 'void *', please use explicit cast [bugprone-multi-level-implicit-pointer-conversion]

PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); ^

github-actions · 2025-10-13T22:17:52Z

src/solver/impls/snes/snes.cxx

  SNESSolver* s;
-  int ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: no header providing "PCShellGetContext" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscpc.h>

github-actions · 2025-10-13T22:17:52Z

src/solver/impls/snes/snes.cxx

  SNESSolver* s;
-  int ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: no header providing "PetscCall" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscerror.h>

github-actions · 2025-10-13T22:17:52Z

src/solver/impls/snes/snes.cxx

                PetscInt col = ind2 + j;
-                ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES);
-                CHKERRQ(ierr);
+                PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES));


warning: no header providing "INSERT_VALUES" is directly included [misc-include-cleaner]

PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES)); ^

github-actions · 2025-10-13T22:17:54Z

src/solver/impls/snes/snes.cxx

+}
+
+// Strategy based on history of residuals
+BoutReal SNESSolver::updatePseudoTimestep_history_based(BoutReal previous_timestep,


warning: method 'updatePseudoTimestep_history_based' can be made const [readability-make-member-function-const]

src/solver/impls/snes/snes.hxx:150:

- BoutReal current_residual); + BoutReal current_residual) const;

src/solver/impls/snes/snes.cxx:1377:

- , + , const

github-actions · 2025-10-13T22:17:54Z

src/solver/impls/snes/snes.cxx

+      return std::max(0.5 * (pseudo_reduction_factor + 1) * previous_timestep,
+                      dt_min_reset);
+    }
+    // Leave unchanged if residual changes are small


warning: do not use 'else' after 'return' [readability-else-after-return]

Suggested change

// Leave unchanged if residual changes are small

;

{

,

;

if (reduction_ratio > 1.2) {

return std::max(0.5 * (pseudo_reduction_factor + 1) * previous_timestep,

dt_min_reset);

}

github-actions · 2025-10-13T22:17:54Z

src/solver/impls/snes/snes.hxx

 BOUT_ENUM_CLASS(BoutSnesEquationForm, pseudo_transient, rearranged_backward_euler,
                backward_euler, direct_newton);

+BOUT_ENUM_CLASS(BoutPTCStrategy,


warning: enum 'BoutPTCStrategy' uses a larger base type ('int', size: 4 bytes) than necessary for its value set, consider using 'std::uint8_t' (1 byte) as the base type to reduce its size [performance-enum-size]

BOUT_ENUM_CLASS(BoutPTCStrategy, ^

github-actions · 2025-10-13T22:17:54Z

src/solver/impls/snes/snes.hxx

 BOUT_ENUM_CLASS(BoutSnesEquationForm, pseudo_transient, rearranged_backward_euler,
                backward_euler, direct_newton);

+BOUT_ENUM_CLASS(BoutPTCStrategy,


warning: function 'BoutPTCStrategyFromString' can be made static or moved into an anonymous namespace to enforce internal linkage [misc-use-internal-linkage]

BOUT_ENUM_CLASS(BoutPTCStrategy, ^

expanded from here

github-actions · 2025-10-13T22:17:54Z

src/solver/impls/snes/snes.hxx

+                                              BoutReal previous_residual,
+                                              BoutReal current_residual);
+  Field3D pseudo_residual; ///< Diagnostic output
+  Field2D pseudo_residual_2d;


warning: no header providing "Field2D" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.hxx:30:

- #include <bout/build_defines.hxx> + #include "bout/field2d.hxx" + #include <bout/build_defines.hxx>

Moving Jacobian calculations into separate functions.

Add manual sections on adaptive timestepping and PTC method. Minor tidying to address some Clang-Tidy comments.

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 55. Check the log or trigger a new build to see more.

github-actions · 2025-10-14T18:13:13Z

src/solver/impls/snes/snes.cxx

+
+    //////////////////////////////////////////////////
+    // Get the local indices by starting at 0
+    Field3D index = globalIndex(0);


warning: no header providing "Field3D" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include "bout/field3d.hxx"

github-actions · 2025-10-14T18:13:14Z

src/solver/impls/snes/snes.cxx

+
+    output_progress.write("Setting Jacobian matrix sizes\n");
+
+    const int n2d = f2d.size();


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

const int n2d = f2d.size(); ^

github-actions · 2025-10-14T18:13:14Z

src/solver/impls/snes/snes.cxx

+    output_progress.write("Setting Jacobian matrix sizes\n");
+
+    const int n2d = f2d.size();
+    const int n3d = f3d.size();


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

const int n3d = f3d.size(); ^

github-actions · 2025-10-14T18:13:14Z

src/solver/impls/snes/snes.cxx

+    const int n3d = f3d.size();
+
+    // Set size of Matrix on each processor to nlocal x nlocal
+    MatCreate(BoutComm::get(), &Jfd);


warning: no header providing "MatCreate" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscmat.h>

github-actions · 2025-10-14T18:13:14Z

src/solver/impls/snes/snes.cxx

+
+    // Set size of Matrix on each processor to nlocal x nlocal
+    MatCreate(BoutComm::get(), &Jfd);
+    MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);


warning: no header providing "MAT_KEEP_NONZERO_PATTERN" is directly included [misc-include-cleaner]

MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); ^

github-actions · 2025-10-14T18:13:16Z

src/solver/impls/snes/snes.cxx

+    MatMPIAIJSetPreallocation(Jfd, 0, d_nnz.data(), 0, o_nnz.data());
+    MatSeqAIJSetPreallocation(Jfd, 0, d_nnz.data());
+    MatSetUp(Jfd);
+    MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);


warning: no header providing "MAT_NEW_NONZERO_ALLOCATION_ERR" is directly included [misc-include-cleaner]

MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); ^

src/solver/impls/snes/snes.cxx

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 55. Check the log or trigger a new build to see more.

github-actions · 2025-10-14T18:13:24Z

src/solver/impls/snes/snes.cxx

+
+    //////////////////////////////////////////////////
+    // Get the local indices by starting at 0
+    Field3D index = globalIndex(0);


warning: no header providing "Field3D" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include "bout/field3d.hxx"

github-actions · 2025-10-14T18:13:24Z

src/solver/impls/snes/snes.cxx

+
+    output_progress.write("Setting Jacobian matrix sizes\n");
+
+    const int n2d = f2d.size();


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

const int n2d = f2d.size(); ^

github-actions · 2025-10-14T18:13:24Z

src/solver/impls/snes/snes.cxx

+    output_progress.write("Setting Jacobian matrix sizes\n");
+
+    const int n2d = f2d.size();
+    const int n3d = f3d.size();


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

const int n3d = f3d.size(); ^

github-actions · 2025-10-14T18:13:24Z

src/solver/impls/snes/snes.cxx

+    const int n3d = f3d.size();
+
+    // Set size of Matrix on each processor to nlocal x nlocal
+    MatCreate(BoutComm::get(), &Jfd);


warning: no header providing "MatCreate" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscmat.h>

github-actions · 2025-10-14T18:13:24Z

src/solver/impls/snes/snes.cxx

+
+    // Set size of Matrix on each processor to nlocal x nlocal
+    MatCreate(BoutComm::get(), &Jfd);
+    MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);


warning: no header providing "MAT_KEEP_NONZERO_PATTERN" is directly included [misc-include-cleaner]

MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); ^

github-actions · 2025-10-14T18:13:26Z

src/solver/impls/snes/snes.cxx

+    MatMPIAIJSetPreallocation(Jfd, 0, d_nnz.data(), 0, o_nnz.data());
+    MatSeqAIJSetPreallocation(Jfd, 0, d_nnz.data());
+    MatSetUp(Jfd);
+    MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);


warning: no header providing "MAT_NEW_NONZERO_ALLOCATION_ERR" is directly included [misc-include-cleaner]

MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); ^

github-actions · 2025-10-14T18:13:26Z

src/solver/impls/snes/snes.cxx

+    // Mark non-zero entries
+
+    output_progress.write("Marking non-zero Jacobian entries\n");
+    PetscScalar val = 1.0;


warning: no header providing "PetscScalar" is directly included [misc-include-cleaner]

PetscScalar val = 1.0; ^

github-actions · 2025-10-14T18:13:26Z

src/solver/impls/snes/snes.cxx

+    // Mark non-zero entries
+
+    output_progress.write("Marking non-zero Jacobian entries\n");
+    PetscScalar val = 1.0;


warning: variable 'val' of type 'PetscScalar' (aka 'double') can be declared 'const' [misc-const-correctness]

Suggested change

PetscScalar val = 1.0;

PetscScalar const val = 1.0;

src/solver/impls/snes/snes.cxx

github-actions

clang-tidy made some suggestions

github-actions · 2025-10-14T22:21:06Z

src/solver/impls/snes/snes.cxx

+        }
+        // 3D fields
+        for (int z = 0; z < mesh->LocalNz; z++) {
+          int ind = ROUND(index(x, y, z));


warning: variable 'ind' of type 'int' can be declared 'const' [misc-const-correctness]

Suggested change

int ind = ROUND(index(x, y, z));

int const ind = ROUND(index(x, y, z));

github-actions · 2025-10-14T22:21:06Z

src/solver/impls/snes/snes.cxx

+
+            // Star pattern
+            for (const auto& [x_off, y_off] : xy_offsets) {
+              int xi = x + x_off;


warning: variable 'xi' of type 'int' can be declared 'const' [misc-const-correctness]

Suggested change

int xi = x + x_off;

int const xi = x + x_off;

github-actions · 2025-10-14T22:21:06Z

src/solver/impls/snes/snes.cxx

+            // Star pattern
+            for (const auto& [x_off, y_off] : xy_offsets) {
+              int xi = x + x_off;
+              int yi = y + y_off;


warning: variable 'yi' of type 'int' can be declared 'const' [misc-const-correctness]

Suggested change

int yi = y + y_off;

int const yi = y + y_off;

github-actions · 2025-10-14T22:21:06Z

src/solver/impls/snes/snes.cxx

+                // 3D fields on this cell
+                for (int j = 0; j < n3d; j++) {
+                  const PetscInt col = ind2 + j;
+                  PetscErrorCode ierr =


warning: variable 'ierr' of type 'PetscErrorCode' (aka 'int') can be declared 'const' [misc-const-correctness]

Suggested change

PetscErrorCode ierr =

PetscErrorCode const ierr =

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+    {
+      // Test if the matrix is symmetric
+      // Values are 0 or 1 so tolerance (1e-5) shouldn't matter
+      PetscBool symmetric;


warning: variable 'symmetric' is not initialized [cppcoreguidelines-init-variables]

PetscBool symmetric; ^

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+      // Values are 0 or 1 so tolerance (1e-5) shouldn't matter
+      PetscBool symmetric;
+      PetscCall(MatIsSymmetric(Jfd, 1e-5, &symmetric));
+      if (!symmetric) {


warning: implicit conversion 'PetscBool' -> 'bool' [readability-implicit-bool-conversion]

Suggested change

if (!symmetric) {

if (symmetric == 0u) {

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+    if ((*options)["force_symmetric_coloring"]
+            .doc("Modifies coloring matrix to force it to be symmetric")
+            .withDefault<bool>(false)) {
+      Mat Jfd_T;


warning: variable 'Jfd_T' is not initialized [cppcoreguidelines-init-variables]

Suggested change

Mat Jfd_T;

Mat Jfd_T = nullptr;

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+  return PETSC_SUCCESS;
+}
+
+PetscErrorCode SNESSolver::FDJpruneJacobian() {


warning: method 'FDJpruneJacobian' can be made static [readability-convert-member-functions-to-static]

src/solver/impls/snes/snes.hxx:104:

- PetscErrorCode FDJpruneJacobian(); ///< Remove small elements from the Jacobian + static PetscErrorCode FDJpruneJacobian(); ///< Remove small elements from the Jacobian

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+}
+
+PetscErrorCode SNESSolver::FDJpruneJacobian() {
+#if PETSC_VERSION_GE(3, 20, 0)


warning: no header providing "PETSC_VERSION_GE" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscversion.h>

github-actions · 2025-10-14T22:21:07Z

src/solver/impls/snes/snes.cxx

+}
+
+// Strategy based on history of residuals
+BoutReal SNESSolver::updatePseudoTimestep_history_based(BoutReal previous_timestep,


warning: method 'updatePseudoTimestep_history_based' can be made const [readability-make-member-function-const]

src/solver/impls/snes/snes.hxx:155:

- BoutReal current_residual); + BoutReal current_residual) const;

src/solver/impls/snes/snes.cxx:1396:

- , + , const

Experimental features, off by default 1. pid_consider_failures (bool, default: false) If enabled, this reduces the magnitude of timestep increases when recent solves have frequently failed. Makes the solver more cautious so may reduce failures but also may not increase timesteps as quickly. Net effect on performance is likely problem dependent. 2. asinh_vars (bool, default: false) If enabled, evolves asinh() of all variables. This might be a good idea when variable magnitudes vary widely, either due to normalisation or strong variations across the mesh. - Behaves like log() for large values - Close to linear for small values, so the solver doesn't over-resolve e.g low density regions - Retains sign, unlike log: asinh(-x) = -asinh(x)

github-actions

clang-tidy made some suggestions

github-actions · 2025-10-15T23:22:12Z

src/solver/impls/snes/snes.cxx

+
+    if (asinh_vars) {
+      // Evolving asinh(vars)
+      PetscInt size;


warning: variable 'size' is not initialized [cppcoreguidelines-init-variables]

Suggested change

PetscInt size;

PetscInt size = 0;

github-actions · 2025-10-15T23:22:12Z

src/solver/impls/snes/snes.cxx

+    }
+
+    if (asinh_vars) {
+      PetscInt size;


warning: variable 'size' is not initialized [cppcoreguidelines-init-variables]

src/solver/impls/snes/snes.cxx:1223:

- ; + = 0;

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

+
+    const BoutReal* xdata = nullptr;
+    PetscCall(VecGetArrayRead(scaled_x, &xdata));
+    load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast to remove const qualifier [cppcoreguidelines-pro-type-const-cast]

; ^

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

+}
+
+// Strategy based on history of residuals
+BoutReal SNESSolver::updatePseudoTimestep_history_based(BoutReal previous_timestep,


warning: method 'updatePseudoTimestep_history_based' can be made const [readability-make-member-function-const]

src/solver/impls/snes/snes.hxx:155:

- BoutReal current_residual); + BoutReal current_residual) const;

src/solver/impls/snes/snes.cxx:1421:

- , + , const

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

  }

+  if (asinh_vars) {
+    PetscInt size;


warning: variable 'size' is not initialized [cppcoreguidelines-init-variables]

src/solver/impls/snes/snes.cxx:1491:

- ; + = 0;

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

+  const BoutReal* xdata = nullptr;
+  PetscCall(VecGetArrayRead(scaled_x, &xdata));
+  // const_cast needed due to load_vars API. Not writing to xdata.
+  load_vars(const_cast<BoutReal*>(xdata));


warning: do not use const_cast to remove const qualifier [cppcoreguidelines-pro-type-const-cast]

. ^

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

+    // du/dt = dvar/dt * du/dvar
+    //
+    // du/var = 1 / sqrt(var^2 + scale^2)
+    PetscInt size;


warning: variable 'size' is not initialized [cppcoreguidelines-init-variables]

src/solver/impls/snes/snes.cxx:1530:

- ; + = 0;

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.cxx

+  if (pid_consider_failures && (fac > 1.0)) {
+    // Reduce aggressiveness if recent steps have failed often
+    fac = pow(fac, std::max(0.3, 1.0 - 2.0 * recent_failure_rate));
+  }


warning: '*' has higher precedence than '-'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]

Suggested change

}

often

ate));()

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.hxx

  BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional)
+  bool pid_consider_failures; ///< Reduce timestep increases if recent solves have failed
+  BoutReal recent_failure_rate;            ///< Rolling average of recent failure rate
+  const BoutReal inv_failure_window = 0.1; ///< 1 / number of recent solves


warning: member 'inv_failure_window' of type 'const BoutReal' (aka 'const double') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

const BoutReal inv_failure_window = 0.1; ///< 1 / number of recent solves ^

github-actions · 2025-10-15T23:22:13Z

src/solver/impls/snes/snes.hxx

  Vec scaled_x;            ///< The values passed to the user RHS
+
+  bool asinh_vars; ///< Evolve asinh(vars) to compress magnitudes while preserving signs
+  const BoutReal asinh_scale = 1e-5; // Scale below which asinh response becomes ~linear


warning: member 'asinh_scale' of type 'const BoutReal' (aka 'const double') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

const BoutReal asinh_scale = 1e-5; // Scale below which asinh response becomes ~linear ^

github-actions bot reviewed Oct 12, 2025

View reviewed changes

bendudson added 2 commits October 13, 2025 10:16

snes: Remove pseudo_time switch, use equation_form=pseudo_transient

e9b8d12

Rather than having a separate switch, just use equation_form to select the Pseudo-Transient Continuation method.

snes: Refactoring, split rhs_function from snes_function

a024d58

The rhs_function returns the time derivatives (residuals) that are used in the PTC method. The snes_function calls rhs_function, then modifies the function depending on equation_form.

github-actions bot reviewed Oct 13, 2025

View reviewed changes

bendudson added 3 commits October 13, 2025 12:23

snes: Add pseudo_strategy option

37187e3

Switches between different strategies for choosing cell time steps.

snes: Replace CHKERRQ with PetscCall macro

07c2882

Minor tidying

snes: Split PTC implementation into separate functions

cc58791

Trying to simplify the `run()` function by splitting out pieces into separate functions.

github-actions bot reviewed Oct 13, 2025

View reviewed changes

bendudson added 2 commits October 13, 2025 19:45

snes: Refactoring Jacobian coloring and pruning

dcbdc5d

Moving Jacobian calculations into separate functions.

snes: Expand manual and minor tidying

48c0b18

Add manual sections on adaptive timestepping and PTC method. Minor tidying to address some Clang-Tidy comments.

github-actions bot reviewed Oct 14, 2025

View reviewed changes

snes: Add missing returns, clang tidying

ca04bc8

github-actions bot reviewed Oct 14, 2025

View reviewed changes

bendudson added 2 commits October 15, 2025 15:39

Merge branch 'next' into snes-pseudo-time

ed642ee

github-actions bot reviewed Oct 15, 2025

View reviewed changes

	#include <bout/petsc_interface.hxx>
	#include <bout/petsc_interface.hxx>


		output_progress.write("Setting Jacobian matrix sizes\n");

		const int n2d = f2d.size();

	int ind = ROUND(index(x, y, z));
	int const ind = ROUND(index(x, y, z));

snes solver: Pseudo Transient Continuation method #3181

Are you sure you want to change the base?

snes solver: Pseudo Transient Continuation method #3181

Uh oh!

Conversation

bendudson commented Oct 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 12, 2025

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions bot left a comment

bendudson commented Oct 12, 2025 •

edited

Loading