Fix indexing issue for output monitoring

include/bout/solver.hxx

@@ -451,19 +451,19 @@ protected:
   ///
   /// There are two important things to note about how \p iter is
   /// passed along to each monitor:
-  /// - The solvers all start their iteration numbering from zero, so the
-  ///   initial state is calculated at \p iter = -1
+  /// - The initial state is written at \p iter = 0, and solver output
+  ///   steps are numbered from 1 to NOUT
   /// - Secondly, \p iter is passed along to each monitor *relative to
   ///   that monitor's period*
   ///
   /// In practice, this means that each monitor is called like:
   ///
   ///     monitor->call(solver, simulation_time,
-  ///                   ((iter + 1) / monitor->period) - 1,
+  ///                   iter / monitor->period,
   ///                   NOUT / monitor->period);
   ///
-  /// e.g. for a monitor with period 10, passing \p iter = 9 will
-  /// result in it being called with a value of `(9 + 1)/10 - 1 == 0`
+  /// e.g. for a monitor with period 10, passing \p iter = 10 will
+  /// result in it being called with a value of `10/10 == 1`
   int call_monitors(BoutReal simtime, int iter, int NOUT);
 
   /// Should timesteps be monitored?

manual/sphinx/user_docs/time_integration.rst

@@ -33,7 +33,7 @@ needed to make the solver available.
 
 .. _tab-solvers:
 .. table:: Available time integration solvers
-	   
+
    +---------------+-----------------------------------------+------------------------+
    | Name          | Description                             | Compile options        |
    +===============+=========================================+========================+
@@ -68,7 +68,7 @@ given in table :numref:`tab-solveropts`.
 
 .. _tab-solveropts:
 .. table:: Time integration solver options
-	   
+
    +--------------------------+--------------------------------------------+-------------------------------------+
    | Option                   | Description                                | Solvers used                        |
    +==========================+============================================+=====================================+
@@ -1278,7 +1278,9 @@ implement the outputMonitor method of PhysicsModel::
     int outputMonitor(BoutReal simtime, int iter, int nout)
 
 The first input is the current simulation time, the second is the output
-number, and the last is the total number of outputs requested.
+number, and the last is the total number of outputs requested. If an initial
+dump is written, it is output number ``0``. Solver output steps are numbered
+from ``1`` to ``nout``, so ``iter == nout`` indicates the final output.
 This method is called by a monitor object PhysicsModel::modelMonitor, which
 writes the restart files at the same time. You can change the frequency at which
 the monitor is called by calling, in PhysicsModel::init::
@@ -1303,7 +1305,9 @@ returns an int::
 
 The first input is the solver object, the second is the current
 simulation time, the third is the output number, and the last is the
-total number of outputs requested. To get the solver to call this
+total number of outputs requested. As for ``outputMonitor()``, output number
+``0`` is reserved for the initial dump when it is written, and solver output
+steps are numbered from ``1`` to ``NOUT``. To get the solver to call this
 function every output time, define a `MyOutputMonitor` object as a member of your
 PhysicsModel::
 

src/solver/impls/adams_bashforth/adams_bashforth.cxx

@@ -355,7 +355,7 @@ int AdamsBashforthSolver::run() {
   [[maybe_unused]] int nwasted = 0;
   [[maybe_unused]] int nwasted_following_fail = 0;
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     bool running = true;

src/solver/impls/arkode/arkode.cxx

@@ -510,7 +510,7 @@ int ArkodeSolver::run() {
     throw BoutException("ArkodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());

src/solver/impls/cvode/cvode.cxx

@@ -489,7 +489,7 @@ int CvodeSolver::run() {
     throw BoutException("CvodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());

src/solver/impls/euler/euler.cxx

@@ -66,7 +66,7 @@ int EulerSolver::init() {
 
 int EulerSolver::run() {
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     bool running = true;

src/solver/impls/ida/ida.cxx

@@ -236,7 +236,7 @@ int IdaSolver::run() {
     throw BoutException("IdaSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());

src/solver/impls/imex-bdf2/imex-bdf2.cxx

@@ -751,7 +751,7 @@ int IMEXBDF2::run() {
 
   int internalCounter = 0; // Cumulative number of successful internal iterations
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal cumulativeTime = 0.;
     int counter = 0; // How many iterations in this output step
 

src/solver/impls/petsc/petsc.cxx

@@ -64,10 +64,10 @@
 
 class ColoringStencil {
 private:
-  bool static isInSquare(int const i, int const j, int const n_square) {
+  bool static isInSquare(const int i, const int j, const int n_square) {
     return std::abs(i) <= n_square && std::abs(j) <= n_square;
   }
-  bool static isInCross(int const i, int const j, int const n_cross) {
+  bool static isInCross(const int i, const int j, const int n_cross) {
     if (i == 0) {
       return std::abs(j) <= n_cross;
     }
@@ -76,7 +76,7 @@ class ColoringStencil {
     }
     return false;
   }
-  bool static isInTaxi(int const i, int const j, int const n_taxi) {
+  bool static isInTaxi(const int i, const int j, const int n_taxi) {
     return std::abs(i) + std::abs(j) <= n_taxi;
   }
 
@@ -234,7 +234,7 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo
     s->load_vars(const_cast<BoutReal*>(x));
     PetscCall(VecRestoreArrayRead(interpolatedX, &x));
 
-    if (s->call_monitors(output_time, i++, s->getNumberOutputSteps()) != 0) {
+    if (s->call_monitors(output_time, ++i, s->getNumberOutputSteps()) != 0) {
       PetscFunctionReturn(1);
     }
 
@@ -614,7 +614,7 @@ int PetscSolver::init() {
         n_taxi = 2;
       }
 
-      auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
+      const auto xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
       {
         // This is ugly but can't think of a better and robust way to
         // count the non-zeros for some arbitrary stencil
@@ -734,7 +734,7 @@ int PetscSolver::init() {
       // Mark non-zero entries
 
       output_progress.write("Marking non-zero Jacobian entries\n");
-      PetscScalar const val = 1.0;
+      const PetscScalar val = 1.0;
       for (int x = mesh->xstart; x <= mesh->xend; x++) {
         for (int y = mesh->ystart; y <= mesh->yend; y++) {
 
@@ -753,21 +753,21 @@ int PetscSolver::init() {
                 continue;
               }
 
-              int const ind2 = ROUND(index(xi, yi, 0));
+              const int ind2 = ROUND(index(xi, yi, 0));
               if (ind2 < 0) {
                 continue; // A boundary point
               }
 
               // Depends on all variables on this cell
               for (int j = 0; j < n2d; j++) {
-                PetscInt const col = ind2 + j;
+                const PetscInt col = ind2 + j;
                 PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES));
               }
             }
           }
           // 3D fields
           for (int z = mesh->zstart; z <= mesh->zend; z++) {
-            int const ind = ROUND(index(x, y, z));
+            const int ind = ROUND(index(x, y, z));
 
             for (int i = 0; i < n3d; i++) {
               PetscInt row = ind + i;
@@ -777,7 +777,7 @@ int PetscSolver::init() {
 
               // Depends on 2D fields
               for (int j = 0; j < n2d; j++) {
-                PetscInt const col = ind0 + j;
+                const PetscInt col = ind0 + j;
                 PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES));
               }
 
@@ -802,7 +802,7 @@ int PetscSolver::init() {
 
                   // 3D fields on this cell
                   for (int j = 0; j < n3d; j++) {
-                    PetscInt const col = ind2 + j;
+                    const PetscInt col = ind2 + j;
                     int ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES);
 
                     if (ierr != 0) {

src/solver/impls/power/power.cxx

@@ -45,7 +45,7 @@ int PowerSolver::run() {
   // Make sure that f0 has a norm of 1
   divide(f0, norm(f0));
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
 
     load_vars(std::begin(f0));
     run_rhs(curtime);

src/solver/impls/pvode/pvode.cxx

@@ -5,7 +5,7 @@
  * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
  *
  * Contact Ben Dudson, bd512@york.ac.uk
- * 
+ *
  * This file is part of BOUT++.
  *
  * BOUT++ is free software: you can redistribute it and/or modify
@@ -20,7 +20,7 @@
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
- * 
+ *
  **************************************************************************/
 
 #include "bout/build_defines.hxx"
@@ -195,8 +195,8 @@ int PvodeSolver::init() {
   BoutReal* udata = N_VDATA(u);
   save_vars(udata);
 
-  /* Call CVodeMalloc to initialize CVODE: 
-     
+  /* Call CVodeMalloc to initialize CVODE:
+
      neq     is the problem size = number of equations
      f       is the user's right hand side function in y'=f(t,y)
      T0      is the initial time
@@ -294,7 +294,7 @@ int PvodeSolver::run() {
     throw BoutException("PvodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());

src/solver/impls/rk3-ssp/rk3-ssp.cxx

@@ -61,7 +61,7 @@ int RK3SSP::init() {
 
 int RK3SSP::run() {
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;

src/solver/impls/rk4/rk4.cxx

@@ -73,7 +73,7 @@ int RK4Solver::init() {
 }
 
 int RK4Solver::run() {
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;

src/solver/impls/rkgeneric/rkgeneric.cxx

@@ -80,7 +80,7 @@ void RKGenericSolver::resetInternalFields() {
 }
 
 int RKGenericSolver::run() {
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;

src/solver/impls/slepc/slepc.cxx

@@ -579,7 +579,7 @@ void SlepcSolver::monitor(PetscInt its, PetscInt nconv, PetscScalar eigr[],
   static bool first = true;
   if (eigenValOnly && first) {
     first = false;
-    resetIterationCounter();
+    resetIterationCounter(1);
   }
 
   // Temporary eigenvalues, converted from the SLEPc eigenvalues
@@ -701,7 +701,7 @@ void SlepcSolver::analyseResults() {
   output << "Converged eigenvalues :\n"
             "\tIndex\tSlepc eig (mag.)\t\t\tBOUT eig (mag.)\n";
 
-  resetIterationCounter();
+  resetIterationCounter(1);
 
   // Declare and create vectors to store eigenfunctions
   Vec vecReal, vecImag;

src/solver/impls/snes/snes.cxx

@@ -15,8 +15,8 @@
 #include <bout/output_bout_types.hxx>
 #include <bout/petsc_interface.hxx>
 #include <bout/solver.hxx>
-#include <bout/utils.hxx>
 #include <bout/unused.hxx>
+#include <bout/utils.hxx>
 
 #include <algorithm>
 #include <cmath>
@@ -35,10 +35,10 @@
 
 class ColoringStencil {
 private:
-  bool static isInSquare(int const i, int const j, int const n_square) {
+  bool static isInSquare(const int i, const int j, const int n_square) {
     return std::abs(i) <= n_square && std::abs(j) <= n_square;
   }
-  bool static isInCross(int const i, int const j, int const n_cross) {
+  bool static isInCross(const int i, const int j, const int n_cross) {
     if (i == 0) {
       return std::abs(j) <= n_cross;
     }
@@ -47,7 +47,7 @@ class ColoringStencil {
     }
     return false;
   }
-  bool static isInTaxi(int const i, int const j, int const n_taxi) {
+  bool static isInTaxi(const int i, const int j, const int n_taxi) {
     return std::abs(i) + std::abs(j) <= n_taxi;
   }
 
@@ -161,7 +161,7 @@ PetscErrorCode SNESSolver::FDJinitialise() {
                       .doc("Extent of stencil (taxi-cab norm)")
                       .withDefault<int>((n_square == 0 && n_cross == 0) ? 2 : 0);
 
-    auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
+    const auto xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
     {
       // This is ugly but can't think of a better and robust way to
       // count the non-zeros for some arbitrary stencil
@@ -543,9 +543,9 @@ SNESSolver::SNESSolver(Options* opts)
       pseudo_alpha((*options)["pseudo_alpha"]
                        .doc("Sets timestep using dt = alpha / residual")
                        .withDefault(100. * atol * timestep)),
-      pseudo_alpha_minimum(
-          (*options)["pseudo_alpha_minimum"].doc("Minimum value of pseudo_alpha")
-          .withDefault(0.1 * pseudo_alpha)),
+      pseudo_alpha_minimum((*options)["pseudo_alpha_minimum"]
+                               .doc("Minimum value of pseudo_alpha")
+                               .withDefault(0.1 * pseudo_alpha)),
       pseudo_growth_factor((*options)["pseudo_growth_factor"]
                                .doc("PTC growth factor on success")
                                .withDefault(1.1)),
@@ -886,7 +886,7 @@ int SNESSolver::run() {
 
   BoutReal target = simtime;
   recent_failure_rate = 0.0;
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     target += getOutputTimestep();
 
     bool looping = true;
@@ -1191,10 +1191,10 @@ int SNESSolver::run() {
 
       if (equation_form == BoutSnesEquationForm::pseudo_transient) {
         // Adjust pseudo_alpha to globally scale timesteps
-        pseudo_alpha = std::max({
-            updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate,
-                                 max_timestep * atol * 100),
-            pseudo_alpha_minimum});
+        pseudo_alpha =
+            std::max({updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate,
+                                           max_timestep * atol * 100),
+                      pseudo_alpha_minimum});
 
         // Adjust local timesteps
         PetscCall(updatePseudoTimestepping());
@@ -1541,10 +1541,10 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() {
 /// rapid changes in timestep.
 BoutReal SNESSolver::updatePseudoTimestep_inverse_residual(BoutReal previous_timestep,
                                                            BoutReal current_residual) {
-  return std::max({
-      std::min({std::max({pseudo_alpha / current_residual, previous_timestep / 1.5}),
-                1.5 * previous_timestep, max_timestep}),
-      dt_min_reset});
+  return std::max(
+      {std::min({std::max({pseudo_alpha / current_residual, previous_timestep / 1.5}),
+                 1.5 * previous_timestep, max_timestep}),
+       dt_min_reset});
 }
 
 // Strategy based on history of residuals