diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 446acefecb..14df410928 100644 --- a/include/bout/solver.hxx +++ b/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? diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 37bf9a6abf..b12e497e02 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/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:: diff --git a/src/solver/impls/adams_bashforth/adams_bashforth.cxx b/src/solver/impls/adams_bashforth/adams_bashforth.cxx index 00bd8f4ff2..0cd39a5b6d 100644 --- a/src/solver/impls/adams_bashforth/adams_bashforth.cxx +++ b/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; diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 804e67a803..11e3561d0c 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/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()); diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index 4f72e65933..2710534a62 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/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()); diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index 9754a68826..599d556c00 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/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; diff --git a/src/solver/impls/ida/ida.cxx b/src/solver/impls/ida/ida.cxx index 84a0bd4abd..68ddf7f90d 100644 --- a/src/solver/impls/ida/ida.cxx +++ b/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()); diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 4ccf27981b..e30ae4c743 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/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 diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 1e1e05596b..8e2ac92d88 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/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(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) { diff --git a/src/solver/impls/power/power.cxx b/src/solver/impls/power/power.cxx index 5dfccfedda..744bdfe46e 100644 --- a/src/solver/impls/power/power.cxx +++ b/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); diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 38901a2a57..ba4dce5849 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/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 . - * + * **************************************************************************/ #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()); diff --git a/src/solver/impls/rk3-ssp/rk3-ssp.cxx b/src/solver/impls/rk3-ssp/rk3-ssp.cxx index 7c00d44b34..26319eaa9a 100644 --- a/src/solver/impls/rk3-ssp/rk3-ssp.cxx +++ b/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; diff --git a/src/solver/impls/rk4/rk4.cxx b/src/solver/impls/rk4/rk4.cxx index 35b3cc2268..bb7c14fab2 100644 --- a/src/solver/impls/rk4/rk4.cxx +++ b/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; diff --git a/src/solver/impls/rkgeneric/rkgeneric.cxx b/src/solver/impls/rkgeneric/rkgeneric.cxx index edf8ed50e3..20f4cc39d6 100644 --- a/src/solver/impls/rkgeneric/rkgeneric.cxx +++ b/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; diff --git a/src/solver/impls/slepc/slepc.cxx b/src/solver/impls/slepc/slepc.cxx index 2c52cb3b38..9887f184eb 100644 --- a/src/solver/impls/slepc/slepc.cxx +++ b/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; diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 5b28ff413c..1dd5e37c84 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -15,8 +15,8 @@ #include #include #include -#include #include +#include #include #include @@ -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((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 diff --git a/src/solver/impls/split-rk/split-rk.cxx b/src/solver/impls/split-rk/split-rk.cxx index ac4d62ebe2..e5fb0f6a8e 100644 --- a/src/solver/impls/split-rk/split-rk.cxx +++ b/src/solver/impls/split-rk/split-rk.cxx @@ -95,7 +95,7 @@ int SplitRK::init() { int SplitRK::run() { - for (int step = 0; step < getNumberOutputSteps(); step++) { + for (int step = 1; step <= getNumberOutputSteps(); step++) { // Take an output step BoutReal target = simtime + getOutputTimestep(); diff --git a/tests/unit/solver/test_fakesolver.hxx b/tests/unit/solver/test_fakesolver.hxx index 3d96b717e3..c8600e4649 100644 --- a/tests/unit/solver/test_fakesolver.hxx +++ b/tests/unit/solver/test_fakesolver.hxx @@ -19,6 +19,13 @@ public: if ((*options)["throw_run"].withDefault(false)) { throw BoutException("Deliberate exception in FakeSolver::run"); } + if ((*options)["call_final_monitor"].withDefault(false)) { + const int ret = + call_monitors(simtime, getNumberOutputSteps(), getNumberOutputSteps()); + if (ret != 0) { + return ret; + } + } return (*options)["fail_run"].withDefault(0); } bool run_called{false}; diff --git a/tests/unit/solver/test_solver.cxx b/tests/unit/solver/test_solver.cxx index e9a612103c..525d52767f 100644 --- a/tests/unit/solver/test_solver.cxx +++ b/tests/unit/solver/test_solver.cxx @@ -971,6 +971,26 @@ TEST_F(SolverTest, BasicSolve) { EXPECT_TRUE(solver.run_called); } +TEST_F(SolverTest, SolveCleansUpMonitorsAtFinalIteration) { + Options options; + options["call_final_monitor"] = true; + FakeSolver solver{&options}; + + StrictMock monitor; + solver.addMonitor(&monitor); + + NiceMock model{}; + solver.setModel(&model); + + Options::cleanup(); + + EXPECT_CALL(monitor, call(_, _, 0, nout)); + EXPECT_CALL(monitor, call(_, _, nout, nout)); + EXPECT_CALL(monitor, cleanup()); + + solver.solve(nout, timestep); +} + TEST_F(SolverTest, GetRunID) { Options options; FakeSolver solver{&options};