From 1c53604a06b600ac455460c111058f6aab9440d9 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 5 Jun 2026 19:15:12 +0100 Subject: [PATCH 01/13] Accumulated work trying to improve scaling for SNES Has involved refactoring Solver to allow loading/saving from/to arbitrary sets of fields. Have also changed how scale_vars works, although I'll probably want to name a different setting for this. What I've made it done doesn't appear to be what it was intended for originally. --- include/bout/solver.hxx | 213 ++++++++++++++++++++++++++++++++- src/solver/impls/snes/snes.cxx | 188 +++++++++++++++++++++-------- src/solver/impls/snes/snes.hxx | 10 ++ src/solver/solver.cxx | 183 +--------------------------- 4 files changed, 365 insertions(+), 229 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 2e0da05224..90b2e5a59c 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -92,7 +92,8 @@ constexpr auto SOLVERIMEXBDF2 = "imexbdf2"; constexpr auto SOLVERSNES = "snes"; constexpr auto SOLVERRKGENERIC = "rkgeneric"; -enum class SOLVER_VAR_OP { LOAD_VARS, LOAD_DERIVS, SET_ID, SAVE_VARS, SAVE_DERIVS }; +enum class FieldCategories { VARS, DERIVS, MMS }; +enum class SOLVER_VAR_OP { LOAD, SET_ID, SAVE }; /// A type to set where in the list monitors are added enum class MonitorPosition { BACK, FRONT }; @@ -391,6 +392,79 @@ protected: std::string description{""}; /// Description of what the variable is }; + /// A structure for iterating over fields + template + struct VarIterator + { + using iterator_category = std::random_access_iterator_tag; + using value_type = T; + using pointer = T*; + using reference = T&; + using underlying_iterator = std::vector>::iterator; + using difference_type = std::iterator_traits::difference_type; + + VarIterator(underlying_iterator _it) : it(std::move(_it)) {} + + reference operator*() const { + switch (C) { + case(FieldCategories::VARS): + return *(it->var); + case(FieldCategories::DERIVS): + return *(it->F_var); + case(FieldCategories::MMS): + return *(it->MMS_err); + } + } + pointer operator->() const { + switch (C) { + case(FieldCategories::VARS): + return it->var; + case(FieldCategories::DERIVS): + return it->F_var; + case(FieldCategories::MMS): + return it->MMS_err.get(); + } + } + + VarIterator& operator++() { it++; return *this; } + VarIterator operator++(int) { VarIterator tmp = *this; ++(*this); return tmp; } + VarIterator& operator+=(difference_type n) { it += n; return *this; } + VarIterator operator+(difference_type n) const { return VarIterator(it + n); } + friend VarIterator operator+(difference_type n, const VarIterator& a) { return VarIterator(n + a.it); } + + VarIterator& operator--() { it--; return *this; } + VarIterator operator--(int) { VarIterator tmp = *this; --(*this); return tmp; } + VarIterator& operator-=(difference_type n) { it -= n; return *this; } + VarIterator operator-(difference_type n) const { return VarIterator(it - n); } + friend VarIterator operator-(difference_type n, const VarIterator& a) { return VarIterator(n - a.it); } + + difference_type operator-(const VarIterator& b) { return it - b.it; } + reference operator[](difference_type n) { return *VarIterator(it[n]); } + + friend bool operator== (const VarIterator& a, const VarIterator& b) { return a.it == b.it; } + friend bool operator!= (const VarIterator& a, const VarIterator& b) { return a.it != b.it; } + friend bool operator< (const VarIterator& a, const VarIterator& b) { return a.it < b.it; } + friend bool operator> (const VarIterator& a, const VarIterator& b) { return a.it > b.it; } + friend bool operator<= (const VarIterator& a, const VarIterator& b) { return a.it <= b.it; } + friend bool operator>= (const VarIterator& a, const VarIterator& b) { return a.it >= b.it; } + + private: + underlying_iterator it; + }; + + template + struct VarRange { + using iterator = VarIterator; + + VarRange(std::vector>& vars) : _begin(vars.begin()), _end(vars.end()) {} + + iterator begin() const { return _begin; } + iterator end() const {return _end; } + private: + iterator _begin, _end; + }; + + /// Does \p var represent field \p name? template friend bool operator==(const VarStr& var, const std::string& name) { @@ -519,6 +593,27 @@ protected: /// Get the currently set output timestep BoutReal getOutputTimestep() const { return output_timestep; } + /// Loop over variables (accessed by iterating the f2d_range and f3d_range + /// arguments) and domain. Used for all data operations for + /// consistency + template + void loop_vars(RangeF2D f2d_range, RangeF3D f3d_range, BoutReal* udata, SOLVER_VAR_OP op) { + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; + + int p = 0; // Counter for location in udata array + + // All boundaries + for (const auto& i2d : mesh->getRegion2D("RGN_BNDRY")) { + loop_vars_op(f2d_range, f3d_range, i2d, udata, p, op, true); + } + + // Bulk of points + for (const auto& i2d : mesh->getRegion2D("RGN_NOBNDRY")) { + loop_vars_op(f2d_range, f3d_range, i2d, udata, p, op, false); + } + } + private: /// Generate a random UUID (version 4) and broadcast it to all processors std::string createRunID() const; @@ -573,9 +668,119 @@ private: /// Should be run after user RHS is called void post_rhs(BoutReal t); - /// Loading data from BOUT++ to/from solver - void loop_vars_op(Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, bool bndry); - void loop_vars(BoutReal* udata, SOLVER_VAR_OP op); + /// Perform an operation at a given Ind2D (jx,jy) location, moving + /// data between BOUT++ and solver. This can be done on arbitrary + /// sets of fields (accessed using f2d_range and f3d_range), provided they + /// are the same number, shape, size, etc. as the fields being + /// solved for. + template + void loop_vars_op(RangeF2D f2d_range, RangeF3D f3d_range, Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, bool bndry) { +/************************************************************************** + * Looping over variables + * + * NOTE: This part is very inefficient, and should be replaced ASAP + * Is the interleaving of variables needed or helpful to the solver? + **************************************************************************/ + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; + + int nz = mesh->LocalNz; + + switch (op) { + case SOLVER_VAR_OP::LOAD: { + /// Load variables from IDA into BOUT++ + + // Loop over 2D variables + auto metadata_it = f2d.begin(); + for (auto& field : f2d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + field[i2d] = udata[p]; + p++; + metadata_it++; + } + + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (auto& field : f3d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + field[field.getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; + p++; + metadata_it++; + } + } + break; + } + case SOLVER_VAR_OP::SET_ID: { + /// Set the type of equation (Differential or Algebraic) + + // Loop over 2D variables + for (const auto& metadata : f2d) { + if (bndry && !metadata.evolve_bndry) { + continue; + } + if (metadata.constraint) { + udata[p] = 0; + } else { + udata[p] = 1; + } + p++; + } + + for (int jz = 0; jz < nz; jz++) { + + // Loop over 3D variables + for (const auto& metadata : f3d) { + if (bndry && !metadata.evolve_bndry) { + continue; + } + if (metadata.constraint) { + udata[p] = 0; + } else { + udata[p] = 1; + } + p++; + } + } + + break; + } + case SOLVER_VAR_OP::SAVE: { + /// Save variables from BOUT++ into IDA (only used at start of simulation) + + // Loop over 2D variables + auto metadata_it = f2d.begin(); + for (const auto& field : f2d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + udata[p] = field[i2d]; + p++; + metadata_it++; + } + + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (const auto& field : f3d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + udata[p] = field[field.getMesh()->ind2Dto3D(i2d, jz)]; + p++; + metadata_it++; + } + } + break; + } + } +} /// Check if a variable has already been added bool varAdded(const std::string& name); diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 1dd5e37c84..d0609b74cf 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -632,6 +632,9 @@ SNESSolver::SNESSolver(Options* opts) scale_vars((*options)["scale_vars"] .doc("Scale variables (Jacobian column scaling)?") .withDefault(false)), + rescale_period((*options)["rescale_period"] + .doc("Number of timesteps before recalculating variabl secaling") + .withDefault(100)), asinh_vars((*options)["asinh_vars"] .doc("Apply asinh() to all variables?") .withDefault(false)) {} @@ -655,6 +658,16 @@ int SNESSolver::init() { output_info.write("\t3d fields = {:d}, 2d fields = {:d} neq={:d}, local_N={:d}\n", n3Dvars(), n2Dvars(), neq, nlocal); + // Initialise fields for storing residual of nonlinear solves + if (diagnose) { + for (const auto& f: f2d) { + resid_2d.emplace_back(emptyFrom(*f.var)); + } + for (const auto& f: f3d) { + resid_3d.emplace_back(emptyFrom(*f.var)); + } + } + // Initialise PETSc components // Vectors @@ -665,6 +678,9 @@ int SNESSolver::init() { PetscCall(VecDuplicate(snes_x, &snes_f)); PetscCall(VecDuplicate(snes_x, &x0)); + if (diagnose) { + PetscCall(VecDuplicate(snes_f, &f0)); + } if ((equation_form == BoutSnesEquationForm::rearranged_backward_euler) || (equation_form == BoutSnesEquationForm::pseudo_transient)) { @@ -755,6 +771,9 @@ int SNESSolver::init() { // Used so that the timestep does not have to be adjusted, // because that would require updating the preconditioner. PetscCall(VecDuplicate(snes_x, &output_x)); + if (diagnose) { + PetscCall(VecDuplicate(snes_f, &output_f)); + } // Initialize the Finite Difference Jacobian PetscCall(FDJinitialise()); @@ -856,6 +875,56 @@ int SNESSolver::init() { return 0; } +PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { + // Individual variable scaling + // Note: If variables are rescaled then the Jacobian columns + // need to be scaled or recalculated + int istart, iend; + VecGetOwnershipRange(snes_x, &istart, &iend); + + // Take ownership of snes_x and var_scaling_factors data + PetscScalar* snes_x_data = nullptr; + PetscCall(VecGetArray(snes_x, &snes_x_data)); + PetscScalar* x1_data; + PetscCall(VecGetArray(x1, &x1_data)); + PetscScalar* var_scaling_factors_data; + PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); + + // Normalise each value in the state + // Limit normalisation so scaling factor is never smaller than rtol + // NOTE: The normalisation value is being progressively updated, hence why we multiply var_scaling_scaling_factors by the new norm + for (int i = 0; i < iend - istart; ++i) { + // THIS LOOKS LIKE IT COULD INTRODUCE A SHARP DISCONTINUITY. COULD THAT BE A SOURCE OF NUMERICAL PROBLEMS? MAYBE, BUT NOTHING TOO SIGNIFICANT THOUGH + const PetscScalar norm = + BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); + snes_x_data[i] /= norm; // FIXME: Is this right? In most cases this means the variable is just 1. Yes, because this is just the guess for the next solution + x1_data[i] /= norm; // Update history for predictor + var_scaling_factors_data[i] *= norm; + } + + // Restore vector underlying data + PetscCall(VecRestoreArray(var_scaling_factors, &var_scaling_factors_data)); + PetscCall(VecRestoreArray(x1, &x1_data)); + PetscCall(VecRestoreArray(snes_x, &snes_x_data)); + + if (diagnose) { + // Print maximum and minimum scaling factors + PetscReal max_scale, min_scale; + VecMax(var_scaling_factors, nullptr, &max_scale); + VecMin(var_scaling_factors, nullptr, &min_scale); + output.write("Var scaling: {} -> {}\n", min_scale, max_scale); + } + + // Force recalculation of the Jacobian + // TODO: Get my head around lagged jacobians so I can check if this is right + // SNESGetLagJacobianPersists(snes, saved_jacobian_persist); + // SNESSetLagJacobianPersists(snes, PETSC_FALSE); + + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + return PETSC_SUCCESS; +} + int SNESSolver::run() { // Set initial guess at the solution from variables { @@ -875,6 +944,11 @@ int SNESSolver::run() { PetscCall(VecRestoreArray(snes_x, &xdata)); } + if (scale_vars) { + int saved_jacobian_lag; + PetscCall(rescale(saved_jacobian_lag)); + } + // Initialise residuals local_residual = 0.0; local_residual_2d = 0.0; @@ -884,6 +958,11 @@ int SNESSolver::run() { output.write("\n Residual: {}\n", global_residual); } + // If saving solver residuals, need to get initial value. + if (diagnose) { + snes_function(snes_x, snes_f, false); + } + BoutReal target = simtime; recent_failure_rate = 0.0; for (int s = 1; s <= getNumberOutputSteps(); s++) { @@ -892,6 +971,7 @@ int SNESSolver::run() { bool looping = true; int snes_failures = 0; // Count SNES convergence failures int saved_jacobian_lag = 0; + PetscBool saved_jacobian_persist = PETSC_FALSE; int loop_count = 0; const BoutReal start_global_residual = global_residual; @@ -902,57 +982,16 @@ int SNESSolver::run() { break; // Could happen if step over multiple outputs } - if (scale_vars) { - // Individual variable scaling - // Note: If variables are rescaled then the Jacobian columns - // need to be scaled or recalculated - - if (loop_count % 100 == 0) { - // Rescale state (snes_x) so that all quantities are around 1 - // If quantities are near zero then RTOL is used - int istart, iend; - VecGetOwnershipRange(snes_x, &istart, &iend); - - // Take ownership of snes_x and var_scaling_factors data - PetscScalar* snes_x_data = nullptr; - PetscCall(VecGetArray(snes_x, &snes_x_data)); - PetscScalar* x1_data; - PetscCall(VecGetArray(x1, &x1_data)); - PetscScalar* var_scaling_factors_data; - PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); - - // Normalise each value in the state - // Limit normalisation so scaling factor is never smaller than rtol - for (int i = 0; i < iend - istart; ++i) { - const PetscScalar norm = - BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); - snes_x_data[i] /= norm; - x1_data[i] /= norm; // Update history for predictor - var_scaling_factors_data[i] *= norm; - } - - // Restore vector underlying data - PetscCall(VecRestoreArray(var_scaling_factors, &var_scaling_factors_data)); - PetscCall(VecRestoreArray(x1, &x1_data)); - PetscCall(VecRestoreArray(snes_x, &snes_x_data)); - - if (diagnose) { - // Print maximum and minimum scaling factors - PetscReal max_scale, min_scale; - VecMax(var_scaling_factors, nullptr, &max_scale); - VecMin(var_scaling_factors, nullptr, &min_scale); - output.write("Var scaling: {} -> {}\n", min_scale, max_scale); - } - - // Force recalculation of the Jacobian - SNESGetLagJacobian(snes, &saved_jacobian_lag); - SNESSetLagJacobian(snes, 1); - } + if (scale_vars and loop_count + 1 % rescale_period == 0) { + PetscCall(rescale(saved_jacobian_lag)); } ++loop_count; // Copy the state (snes_x) into initial values (x0) VecCopy(snes_x, x0); + if (diagnose) { + VecCopy(snes_f, f0); + } if (equation_form == BoutSnesEquationForm::pseudo_transient) { // Pseudo-Transient Continuation @@ -991,6 +1030,8 @@ int SNESSolver::run() { // Note: When the timestep is changed the preconditioner needs to be updated // => Step over the output time and interpolate if not matrix free + // FIXME: Is this the most appropriate thing to do when we're taking such massive timesteps? + if (matrix_free) { // Ensure that the timestep goes to the next output time and then stops. // This avoids the need to interpolate @@ -1077,6 +1118,9 @@ int SNESSolver::run() { } // Restore state VecCopy(x0, snes_x); + if (diagnose) { + VecCopy(f0, snes_f); + } // Recalculate the Jacobian if (jacobian_pruned and (snes_failures > 2) and (4 * lin_its > 3 * maxl)) { @@ -1088,6 +1132,7 @@ int SNESSolver::run() { PetscCall(FDJrestoreFromPruning()); } + // **** FIXME: Replace with booleans if (saved_jacobian_lag == 0) { // This triggers a Jacobian recalculation SNESGetLagJacobian(snes, &saved_jacobian_lag); @@ -1108,6 +1153,7 @@ int SNESSolver::run() { continue; // Try again } + // **** FIXME: Replace with booleans if (saved_jacobian_lag != 0) { // Following successful step, reset Jacobian lag // to previous value @@ -1232,11 +1278,19 @@ int SNESSolver::run() { // output_x <- alpha * x0 + (1 - alpha) * output_x VecAXPBY(output_x, alpha, 1. - alpha, x0); - output_time = target; + if (diagnose) { + VecCopy(snes_f, output_f); + VecAXPBY(output_f, alpha, 1 - alpha, f0); + } + + output_time = target; } else { // Timestep was adjusted to hit target output time output_x = snes_x; + if (diagnose) { + output_f = snes_f; + } } // Put the result into variables @@ -1244,6 +1298,7 @@ int SNESSolver::run() { // scaled_x <- output_x * var_scaling_factors PetscCall(VecPointwiseMult(scaled_x, output_x, var_scaling_factors)); } else if (asinh_vars) { + // FIXME: scale_vars and asinh_vars aren't mutually exclusive, are they? PetscCall(VecCopy(output_x, scaled_x)); } else { scaled_x = output_x; @@ -1275,6 +1330,13 @@ int SNESSolver::run() { BoutComm::abort(1); } + if (diagnose) { + BoutReal* resid_data; + PetscCall(VecGetArray(output_f, &resid_data)); + loop_vars(resid_2d, resid_3d, resid_data, SOLVER_VAR_OP::LOAD); + PetscCall(VecRestoreArray(output_f, &resid_data)); + } + if (call_monitors(output_time, s, getNumberOutputSteps()) != 0) { break; // User signalled to quit } @@ -1672,6 +1734,10 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) { } PetscCall(VecRestoreArray(f, &fdata)); + + if (scale_vars) { + PetscCall(VecPointwiseDivide(f, f, var_scaling_factors)); + } return PETSC_SUCCESS; } @@ -1692,6 +1758,9 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // f = (x0 - x)/Δt + f // First calculate x - x0 to minimise floating point issues VecWAXPY(delta_x, -1.0, x0, x); // delta_x = x - x0 + // if (scale_vars) { + // VecPointwiseMult(delta_x, delta_x, var_scaling_factors); + // } VecAXPY(f, -1. / dt, delta_x); // f <- f - delta_x / dt break; } @@ -1700,6 +1769,9 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // except that Δt is a vector // f = (x0 - x)/Δt + f VecWAXPY(delta_x, -1.0, x0, x); + // if (scale_vars) { + // VecPointwiseMult(delta_x, delta_x, var_scaling_factors); + // } VecPointwiseDivide(delta_x, delta_x, dt_vec); // delta_x /= dt VecAXPY(f, -1., delta_x); // f <- f - delta_x break; @@ -1707,6 +1779,10 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { case BoutSnesEquationForm::backward_euler: { // Backward Euler // Set f = x - x0 - Δt*f + // if (scale_vars) { + // VecPointwiseMult(x, x, var_scaling_factors); + // VecPointwiseMult(x0, x0, var_scaling_factors); + // } VecAYPX(f, -dt, x); // f <- x - Δt*f VecAXPY(f, -1.0, x0); // f <- f - x0 break; @@ -1933,6 +2009,24 @@ void SNESSolver::outputVars(Options& output_options, bool save_repeat) { output_options["snes_pseudo_timestep"].assignRepeat(pseudo_timestep, "t", save_repeat, "SNESSolver"); } + + if (diagnose) { + auto meta_2d_it = f2d.begin(); + for (const auto& f : resid_2d) { + const std::string name = "resid_" + meta_2d_it->name; + output_options[name].assignRepeat(f, "t", save_repeat, "SNES"); + output_options[name].attributes["description"] = meta_2d_it->description; + ++meta_2d_it; + } + + auto meta_3d_it = f3d.begin(); + for (const auto& f : resid_3d) { + const std::string name = "resid_" + meta_3d_it->name; + output_options[name].assignRepeat(f, "t", save_repeat, "SNES"); + output_options[name].attributes["description"] = meta_3d_it->description; + ++meta_3d_it; + } + } } #endif // BOUT_HAS_PETSC diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index fbcf3baddf..2b5e0eb30d 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -115,6 +115,10 @@ private: PetscErrorCode FDJpruneJacobian(); ///< Remove small elements from the Jacobian PetscErrorCode FDJrestoreFromPruning(); ///< Restore Jacobian to original pattern + /// Rescale state (snes_x) so that all quantities are around 1. If + /// quantities are near zero then RTOL is used. + PetscErrorCode rescale(int& saved_jacobian_lag); + /// Call the physics model RHS function /// /// @param[in] x The state vector. Will be scaled if scale_vars=true @@ -220,8 +224,10 @@ private: Vec snes_f; ///< Used by SNES to store function Vec snes_x; ///< Result of SNES Vec x0; ///< Solution at start of current timestep + Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) Vec delta_x; ///< Change in solution Vec output_x; ///< Solution to output. Used if interpolating. + Vec output_f; ///< Residual to output, if diagnose == true. Used if interpolating. bool predictor; ///< Use linear predictor? Vec x1; ///< Previous solution @@ -260,11 +266,15 @@ private: Vec jac_row_inv_norms; ///< 1 / Norm of the rows of the Jacobian bool scale_vars; ///< Scale individual variables? + int rescale_period; ///< How many time-steps before rescaling variables Vec var_scaling_factors; ///< Factors to multiply variables when passing to user 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 + + std::vector resid_2d; ///< Storage for residuals of SNES solve, unpacked from snes_f + std::vector resid_3d; ///< Storage for residuals of SNES solve, unpacked from snes_f }; #else diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index e105043c7c..0afb21f528 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1028,179 +1028,6 @@ std::unique_ptr Solver::create(const SolverType& type, Options* opts) { return SolverFactory::getInstance().create(type, opts); } -/************************************************************************** - * Looping over variables - * - * NOTE: This part is very inefficient, and should be replaced ASAP - * Is the interleaving of variables needed or helpful to the solver? - **************************************************************************/ - -/// Perform an operation at a given Ind2D (jx,jy) location, moving data between BOUT++ and CVODE -void Solver::loop_vars_op(Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, - bool bndry) { - // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; - - int nz = mesh->LocalNz; - - switch (op) { - case SOLVER_VAR_OP::LOAD_VARS: { - /// Load variables from IDA into BOUT++ - - // Loop over 2D variables - for (const auto& f : f2d) { - if (bndry && !f.evolve_bndry) { - continue; - } - (*f.var)[i2d] = udata[p]; - p++; - } - - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& f : f3d) { - if (bndry && !f.evolve_bndry) { - continue; - } - (*f.var)[f.var->getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; - p++; - } - } - break; - } - case SOLVER_VAR_OP::LOAD_DERIVS: { - /// Load derivatives from IDA into BOUT++ - /// Used for preconditioner - - // Loop over 2D variables - for (const auto& f : f2d) { - if (bndry && !f.evolve_bndry) { - continue; - } - (*f.F_var)[i2d] = udata[p]; - p++; - } - - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& f : f3d) { - if (bndry && !f.evolve_bndry) { - continue; - } - (*f.F_var)[f.F_var->getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; - p++; - } - } - - break; - } - case SOLVER_VAR_OP::SET_ID: { - /// Set the type of equation (Differential or Algebraic) - - // Loop over 2D variables - for (const auto& f : f2d) { - if (bndry && !f.evolve_bndry) { - continue; - } - if (f.constraint) { - udata[p] = 0; - } else { - udata[p] = 1; - } - p++; - } - - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& f : f3d) { - if (bndry && !f.evolve_bndry) { - continue; - } - if (f.constraint) { - udata[p] = 0; - } else { - udata[p] = 1; - } - p++; - } - } - - break; - } - case SOLVER_VAR_OP::SAVE_VARS: { - /// Save variables from BOUT++ into IDA (only used at start of simulation) - - // Loop over 2D variables - for (const auto& f : f2d) { - if (bndry && !f.evolve_bndry) { - continue; - } - udata[p] = (*f.var)[i2d]; - p++; - } - - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& f : f3d) { - if (bndry && !f.evolve_bndry) { - continue; - } - udata[p] = (*f.var)[f.var->getMesh()->ind2Dto3D(i2d, jz)]; - p++; - } - } - break; - } - /// Save time-derivatives from BOUT++ into CVODE (returning RHS result) - case SOLVER_VAR_OP::SAVE_DERIVS: { - - // Loop over 2D variables - for (const auto& f : f2d) { - if (bndry && !f.evolve_bndry) { - continue; - } - udata[p] = (*f.F_var)[i2d]; - p++; - } - - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& f : f3d) { - if (bndry && !f.evolve_bndry) { - continue; - } - udata[p] = (*f.F_var)[f.F_var->getMesh()->ind2Dto3D(i2d, jz)]; - p++; - } - } - break; - } - } -} - -/// Loop over variables and domain. Used for all data operations for consistency -void Solver::loop_vars(BoutReal* udata, SOLVER_VAR_OP op) { - // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; - - int p = 0; // Counter for location in udata array - - // All boundaries - for (const auto& i2d : mesh->getRegion2D("RGN_BNDRY")) { - loop_vars_op(i2d, udata, p, op, true); - } - - // Bulk of points - for (const auto& i2d : mesh->getRegion2D("RGN_NOBNDRY")) { - loop_vars_op(i2d, udata, p, op, false); - } -} - void Solver::load_vars(BoutReal* udata) { // Make sure data is allocated for (const auto& f : f2d) { @@ -1211,7 +1038,7 @@ void Solver::load_vars(BoutReal* udata) { f.var->setLocation(f.location); } - loop_vars(udata, SOLVER_VAR_OP::LOAD_VARS); + loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); // Mark each vector as either co- or contra-variant @@ -1233,7 +1060,7 @@ void Solver::load_derivs(BoutReal* udata) { f.F_var->setLocation(f.location); } - loop_vars(udata, SOLVER_VAR_OP::LOAD_DERIVS); + loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); // Mark each vector as either co- or contra-variant @@ -1275,7 +1102,7 @@ void Solver::save_vars(BoutReal* udata) { } } - loop_vars(udata, SOLVER_VAR_OP::SAVE_VARS); + loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::SAVE); } void Solver::save_derivs(BoutReal* dudata) { @@ -1305,10 +1132,10 @@ void Solver::save_derivs(BoutReal* dudata) { } } - loop_vars(dudata, SOLVER_VAR_OP::SAVE_DERIVS); + loop_vars(VarRange(f2d), VarRange(f3d), dudata, SOLVER_VAR_OP::SAVE); } -void Solver::set_id(BoutReal* udata) { loop_vars(udata, SOLVER_VAR_OP::SET_ID); } +void Solver::set_id(BoutReal* udata) { loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::SET_ID); } Field3D Solver::globalIndex(int localStart) { // Use global mesh: FIX THIS! From 769113792c4dd3ca047a1ba3cb77354f89587ca1 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Tue, 16 Jun 2026 15:25:58 +0100 Subject: [PATCH 02/13] Fix rescale period --- 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 d0609b74cf..21b6ce2029 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -982,7 +982,7 @@ int SNESSolver::run() { break; // Could happen if step over multiple outputs } - if (scale_vars and loop_count + 1 % rescale_period == 0) { + if (scale_vars and (loop_count + 1) % rescale_period == 0) { PetscCall(rescale(saved_jacobian_lag)); } ++loop_count; From 13143335ab9c09d0f97ba52a2047f5efae67d4f1 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Wed, 17 Jun 2026 17:30:37 +0100 Subject: [PATCH 03/13] Make sure not to overwrite function eval with time deriv --- src/solver/impls/snes/snes.cxx | 44 +++++++++++++++++++++++++++++----- src/solver/impls/snes/snes.hxx | 1 + 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 21b6ce2029..8b92712a0e 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -680,6 +680,7 @@ int SNESSolver::init() { PetscCall(VecDuplicate(snes_x, &x0)); if (diagnose) { PetscCall(VecDuplicate(snes_f, &f0)); + PetscCall(VecDuplicate(snes_f, &deriv)); } if ((equation_form == BoutSnesEquationForm::rearranged_backward_euler) @@ -982,7 +983,7 @@ int SNESSolver::run() { break; // Could happen if step over multiple outputs } - if (scale_vars and (loop_count + 1) % rescale_period == 0) { + if (scale_vars and loop_count % rescale_period == 0 and (s > 0 or loop_count > 0)) { PetscCall(rescale(saved_jacobian_lag)); } ++loop_count; @@ -1118,6 +1119,7 @@ int SNESSolver::run() { } // Restore state VecCopy(x0, snes_x); + // FIXME: Not sure I need this if (diagnose) { VecCopy(f0, snes_f); } @@ -1198,20 +1200,37 @@ int SNESSolver::run() { } // Copy derivatives back - { + if (diagnose) { + BoutReal* fdata = nullptr; + PetscCall(VecGetArray(deriv, &fdata)); + save_derivs(fdata); + PetscCall(VecRestoreArray(deriv, &fdata)); + // Forward Euler + VecAXPY(snes_x, dt, deriv); + } + else { BoutReal* fdata = nullptr; PetscCall(VecGetArray(snes_f, &fdata)); save_derivs(fdata); PetscCall(VecRestoreArray(snes_f, &fdata)); + // Forward Euler + VecAXPY(snes_x, dt, snes_f); } - // Forward Euler - VecAXPY(snes_x, dt, snes_f); } + // int i; + // PetscReal norm, minval, maxval; + // PetscCall(VecNorm(snes_f, NORM_2, &norm)); + // PetscCall(VecMin(snes_f, &i, &minval)); + // PetscCall(VecMax(snes_f, &i, &maxval)); + // output.write("Max, min, and norm of residual for SNES function: {}, {}, {}\n", maxval, minval, norm); + + simtime += dt; // Update local and global residuals + /// FIXME: This overwrites snes_f with the time derivs! PetscCall(updateResiduals(snes_x)); if (diagnose) { @@ -1280,6 +1299,7 @@ int SNESSolver::run() { VecAXPBY(output_x, alpha, 1. - alpha, x0); if (diagnose) { + // FIXME: Is this a good approach to take here? It's not the *actual* residual of a solve this way. VecCopy(snes_f, output_f); VecAXPBY(output_f, alpha, 1 - alpha, f0); } @@ -1392,12 +1412,20 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { local_residual_2d_prev = copy(local_residual_2d); global_residual_prev = global_residual; + const BoutReal* current_residual = nullptr; + if (diagnose) { + // Call RHS function to get time derivatives + PetscCall(rhs_function(x, deriv, false)); + + // Reading the residual vectors + PetscCall(VecGetArrayRead(deriv, ¤t_residual)); + } else { // Call RHS function to get time derivatives PetscCall(rhs_function(x, snes_f, false)); // Reading the residual vectors - const BoutReal* current_residual = nullptr; PetscCall(VecGetArrayRead(snes_f, ¤t_residual)); + } // Note: The ordering of quantities in the PETSc vectors // depends on the Solver::loop_vars function @@ -1465,8 +1493,12 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { } } + if (diagnose) { // Restore Vec data arrays - PetscCall(VecRestoreArrayRead(snes_f, ¤t_residual)); + PetscCall(VecRestoreArrayRead(deriv, ¤t_residual)); + } else { + PetscCall(VecRestoreArrayRead(snes_f, ¤t_residual)); + } // Global residual metric (RMS) global_residual = std::sqrt(mean(SQ(local_residual), true)); diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 2b5e0eb30d..43e1783b07 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -222,6 +222,7 @@ private: PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function + Vec deriv; ///< Time derivative; only used if diagnose = true, otherwise will store in snes_f Vec snes_x; ///< Result of SNES Vec x0; ///< Solution at start of current timestep Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) From 227bba89301073a8d26a5065311e7bc3f395af34 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 19 Jun 2026 16:30:15 +0100 Subject: [PATCH 04/13] Add dynamic algorithm for rescaling variables in SNES --- src/solver/impls/snes/snes.cxx | 18 +++++++++++++----- src/solver/impls/snes/snes.hxx | 1 + 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 8b92712a0e..d36b0c13e6 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -633,8 +633,11 @@ SNESSolver::SNESSolver(Options* opts) .doc("Scale variables (Jacobian column scaling)?") .withDefault(false)), rescale_period((*options)["rescale_period"] - .doc("Number of timesteps before recalculating variabl secaling") - .withDefault(100)), + .doc("Number of iterations before recalculating variable scaling") + .withDefault(30)), + rescale_threshold((*options)["rescale_threshold"] + .doc("How much change their should be in the norm of the state, before rescaling") + .withDefault(100.)), asinh_vars((*options)["asinh_vars"] .doc("Apply asinh() to all variables?") .withDefault(false)) {} @@ -965,7 +968,9 @@ int SNESSolver::run() { } BoutReal target = simtime; + BoutReal accumulated_change = 0.; recent_failure_rate = 0.0; + int loop_count = 0; for (int s = 1; s <= getNumberOutputSteps(); s++) { target += getOutputTimestep(); @@ -973,7 +978,6 @@ int SNESSolver::run() { int snes_failures = 0; // Count SNES convergence failures int saved_jacobian_lag = 0; PetscBool saved_jacobian_persist = PETSC_FALSE; - int loop_count = 0; const BoutReal start_global_residual = global_residual; do { @@ -983,10 +987,12 @@ int SNESSolver::run() { break; // Could happen if step over multiple outputs } - if (scale_vars and loop_count % rescale_period == 0 and (s > 0 or loop_count > 0)) { + //output.write("Accumulated change to state is: {}", accumulated_change); + if (scale_vars and (accumulated_change > rescale_threshold or loop_count == rescale_period)) { PetscCall(rescale(saved_jacobian_lag)); + accumulated_change = 0.; + loop_count = 0; } - ++loop_count; // Copy the state (snes_x) into initial values (x0) VecCopy(snes_x, x0); @@ -1232,6 +1238,8 @@ int SNESSolver::run() { // Update local and global residuals /// FIXME: This overwrites snes_f with the time derivs! PetscCall(updateResiduals(snes_x)); + accumulated_change += global_residual * dt; + ++loop_count; if (diagnose) { // Gather and print diagnostic information diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 43e1783b07..82f070805a 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -268,6 +268,7 @@ private: bool scale_vars; ///< Scale individual variables? int rescale_period; ///< How many time-steps before rescaling variables + BoutReal rescale_threshold; //< How much change in the state there should be before rescaling Vec var_scaling_factors; ///< Factors to multiply variables when passing to user Vec scaled_x; ///< The values passed to the user RHS From 9fa91a0c4528bde68f78a5a568e40fa85cca5bcc Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 19 Jun 2026 17:04:47 +0100 Subject: [PATCH 05/13] Tidying up (formatting, removing notes I left, etc.) for PR --- include/bout/solver.hxx | 229 +++++++++++++++++++-------------- src/solver/impls/snes/snes.cxx | 164 +++++++++++------------ src/solver/impls/snes/snes.hxx | 17 ++- src/solver/solver.cxx | 17 ++- 4 files changed, 231 insertions(+), 196 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 90b2e5a59c..fa67ac1ca7 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -394,8 +394,7 @@ protected: /// A structure for iterating over fields template - struct VarIterator - { + struct VarIterator { using iterator_category = std::random_access_iterator_tag; using value_type = T; using pointer = T*; @@ -407,46 +406,82 @@ protected: reference operator*() const { switch (C) { - case(FieldCategories::VARS): + case (FieldCategories::VARS): return *(it->var); - case(FieldCategories::DERIVS): + case (FieldCategories::DERIVS): return *(it->F_var); - case(FieldCategories::MMS): + case (FieldCategories::MMS): return *(it->MMS_err); } } pointer operator->() const { switch (C) { - case(FieldCategories::VARS): + case (FieldCategories::VARS): return it->var; - case(FieldCategories::DERIVS): + case (FieldCategories::DERIVS): return it->F_var; - case(FieldCategories::MMS): + case (FieldCategories::MMS): return it->MMS_err.get(); } } - VarIterator& operator++() { it++; return *this; } - VarIterator operator++(int) { VarIterator tmp = *this; ++(*this); return tmp; } - VarIterator& operator+=(difference_type n) { it += n; return *this; } + VarIterator& operator++() { + it++; + return *this; + } + VarIterator operator++(int) { + VarIterator tmp = *this; + ++(*this); + return tmp; + } + VarIterator& operator+=(difference_type n) { + it += n; + return *this; + } VarIterator operator+(difference_type n) const { return VarIterator(it + n); } - friend VarIterator operator+(difference_type n, const VarIterator& a) { return VarIterator(n + a.it); } + friend VarIterator operator+(difference_type n, const VarIterator& a) { + return VarIterator(n + a.it); + } - VarIterator& operator--() { it--; return *this; } - VarIterator operator--(int) { VarIterator tmp = *this; --(*this); return tmp; } - VarIterator& operator-=(difference_type n) { it -= n; return *this; } + VarIterator& operator--() { + it--; + return *this; + } + VarIterator operator--(int) { + VarIterator tmp = *this; + --(*this); + return tmp; + } + VarIterator& operator-=(difference_type n) { + it -= n; + return *this; + } VarIterator operator-(difference_type n) const { return VarIterator(it - n); } - friend VarIterator operator-(difference_type n, const VarIterator& a) { return VarIterator(n - a.it); } + friend VarIterator operator-(difference_type n, const VarIterator& a) { + return VarIterator(n - a.it); + } difference_type operator-(const VarIterator& b) { return it - b.it; } reference operator[](difference_type n) { return *VarIterator(it[n]); } - - friend bool operator== (const VarIterator& a, const VarIterator& b) { return a.it == b.it; } - friend bool operator!= (const VarIterator& a, const VarIterator& b) { return a.it != b.it; } - friend bool operator< (const VarIterator& a, const VarIterator& b) { return a.it < b.it; } - friend bool operator> (const VarIterator& a, const VarIterator& b) { return a.it > b.it; } - friend bool operator<= (const VarIterator& a, const VarIterator& b) { return a.it <= b.it; } - friend bool operator>= (const VarIterator& a, const VarIterator& b) { return a.it >= b.it; } + + friend bool operator==(const VarIterator& a, const VarIterator& b) { + return a.it == b.it; + } + friend bool operator!=(const VarIterator& a, const VarIterator& b) { + return a.it != b.it; + } + friend bool operator<(const VarIterator& a, const VarIterator& b) { + return a.it < b.it; + } + friend bool operator>(const VarIterator& a, const VarIterator& b) { + return a.it > b.it; + } + friend bool operator<=(const VarIterator& a, const VarIterator& b) { + return a.it <= b.it; + } + friend bool operator>=(const VarIterator& a, const VarIterator& b) { + return a.it >= b.it; + } private: underlying_iterator it; @@ -459,12 +494,12 @@ protected: VarRange(std::vector>& vars) : _begin(vars.begin()), _end(vars.end()) {} iterator begin() const { return _begin; } - iterator end() const {return _end; } + iterator end() const { return _end; } + private: iterator _begin, _end; }; - /// Does \p var represent field \p name? template friend bool operator==(const VarStr& var, const std::string& name) { @@ -596,8 +631,9 @@ protected: /// Loop over variables (accessed by iterating the f2d_range and f3d_range /// arguments) and domain. Used for all data operations for /// consistency - template - void loop_vars(RangeF2D f2d_range, RangeF3D f3d_range, BoutReal* udata, SOLVER_VAR_OP op) { + template + void loop_vars(RangeF2D f2d_range, RangeF3D f3d_range, BoutReal* udata, + SOLVER_VAR_OP op) { // Use global mesh: FIX THIS! Mesh* mesh = bout::globals::mesh; @@ -613,7 +649,7 @@ protected: loop_vars_op(f2d_range, f3d_range, i2d, udata, p, op, false); } } - + private: /// Generate a random UUID (version 4) and broadcast it to all processors std::string createRunID() const; @@ -673,69 +709,55 @@ private: /// sets of fields (accessed using f2d_range and f3d_range), provided they /// are the same number, shape, size, etc. as the fields being /// solved for. - template - void loop_vars_op(RangeF2D f2d_range, RangeF3D f3d_range, Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, bool bndry) { -/************************************************************************** + template + void loop_vars_op(RangeF2D f2d_range, RangeF3D f3d_range, Ind2D i2d, BoutReal* udata, + int& p, SOLVER_VAR_OP op, bool bndry) { + /************************************************************************** * Looping over variables * * NOTE: This part is very inefficient, and should be replaced ASAP * Is the interleaving of variables needed or helpful to the solver? **************************************************************************/ - // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; - - int nz = mesh->LocalNz; - - switch (op) { - case SOLVER_VAR_OP::LOAD: { - /// Load variables from IDA into BOUT++ + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; - // Loop over 2D variables - auto metadata_it = f2d.begin(); - for (auto& field : f2d_range) { - if (bndry && !metadata_it->evolve_bndry) { - continue; - } - field[i2d] = udata[p]; - p++; - metadata_it++; - } + int nz = mesh->LocalNz; - for (int jz = 0; jz < nz; jz++) { + switch (op) { + case SOLVER_VAR_OP::LOAD: { + /// Load variables from IDA into BOUT++ - auto metadata_it = f3d.begin(); - // Loop over 3D variables - for (auto& field : f3d_range) { + // Loop over 2D variables + auto metadata_it = f2d.begin(); + for (auto& field : f2d_range) { if (bndry && !metadata_it->evolve_bndry) { continue; } - field[field.getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; + field[i2d] = udata[p]; p++; metadata_it++; } - } - break; - } - case SOLVER_VAR_OP::SET_ID: { - /// Set the type of equation (Differential or Algebraic) - // Loop over 2D variables - for (const auto& metadata : f2d) { - if (bndry && !metadata.evolve_bndry) { - continue; - } - if (metadata.constraint) { - udata[p] = 0; - } else { - udata[p] = 1; + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (auto& field : f3d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + field[field.getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; + p++; + metadata_it++; + } } - p++; + break; } + case SOLVER_VAR_OP::SET_ID: { + /// Set the type of equation (Differential or Algebraic) - for (int jz = 0; jz < nz; jz++) { - - // Loop over 3D variables - for (const auto& metadata : f3d) { + // Loop over 2D variables + for (const auto& metadata : f2d) { if (bndry && !metadata.evolve_bndry) { continue; } @@ -746,41 +768,56 @@ private: } p++; } - } - break; - } - case SOLVER_VAR_OP::SAVE: { - /// Save variables from BOUT++ into IDA (only used at start of simulation) - - // Loop over 2D variables - auto metadata_it = f2d.begin(); - for (const auto& field : f2d_range) { - if (bndry && !metadata_it->evolve_bndry) { - continue; + for (int jz = 0; jz < nz; jz++) { + + // Loop over 3D variables + for (const auto& metadata : f3d) { + if (bndry && !metadata.evolve_bndry) { + continue; + } + if (metadata.constraint) { + udata[p] = 0; + } else { + udata[p] = 1; + } + p++; + } } - udata[p] = field[i2d]; - p++; - metadata_it++; - } - for (int jz = 0; jz < nz; jz++) { + break; + } + case SOLVER_VAR_OP::SAVE: { + /// Save variables from BOUT++ into IDA (only used at start of simulation) - auto metadata_it = f3d.begin(); - // Loop over 3D variables - for (const auto& field : f3d_range) { + // Loop over 2D variables + auto metadata_it = f2d.begin(); + for (const auto& field : f2d_range) { if (bndry && !metadata_it->evolve_bndry) { continue; } - udata[p] = field[field.getMesh()->ind2Dto3D(i2d, jz)]; + udata[p] = field[i2d]; p++; metadata_it++; } + + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (const auto& field : f3d_range) { + if (bndry && !metadata_it->evolve_bndry) { + continue; + } + udata[p] = field[field.getMesh()->ind2Dto3D(i2d, jz)]; + p++; + metadata_it++; + } + } + break; + } } - break; - } } -} /// Check if a variable has already been added bool varAdded(const std::string& name); diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index d36b0c13e6..21757cda90 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -632,12 +632,14 @@ SNESSolver::SNESSolver(Options* opts) scale_vars((*options)["scale_vars"] .doc("Scale variables (Jacobian column scaling)?") .withDefault(false)), - rescale_period((*options)["rescale_period"] - .doc("Number of iterations before recalculating variable scaling") - .withDefault(30)), + rescale_period( + (*options)["rescale_period"] + .doc("Number of iterations before recalculating variable scaling") + .withDefault(30)), rescale_threshold((*options)["rescale_threshold"] - .doc("How much change their should be in the norm of the state, before rescaling") - .withDefault(100.)), + .doc("How much change their should be in the norm of the " + "state, before rescaling") + .withDefault(100.)), asinh_vars((*options)["asinh_vars"] .doc("Apply asinh() to all variables?") .withDefault(false)) {} @@ -663,10 +665,10 @@ int SNESSolver::init() { // Initialise fields for storing residual of nonlinear solves if (diagnose) { - for (const auto& f: f2d) { + for (const auto& f : f2d) { resid_2d.emplace_back(emptyFrom(*f.var)); } - for (const auto& f: f3d) { + for (const auto& f : f3d) { resid_3d.emplace_back(emptyFrom(*f.var)); } } @@ -880,53 +882,48 @@ int SNESSolver::init() { } PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { - // Individual variable scaling - // Note: If variables are rescaled then the Jacobian columns - // need to be scaled or recalculated - int istart, iend; - VecGetOwnershipRange(snes_x, &istart, &iend); - - // Take ownership of snes_x and var_scaling_factors data - PetscScalar* snes_x_data = nullptr; - PetscCall(VecGetArray(snes_x, &snes_x_data)); - PetscScalar* x1_data; - PetscCall(VecGetArray(x1, &x1_data)); - PetscScalar* var_scaling_factors_data; - PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); - - // Normalise each value in the state - // Limit normalisation so scaling factor is never smaller than rtol - // NOTE: The normalisation value is being progressively updated, hence why we multiply var_scaling_scaling_factors by the new norm - for (int i = 0; i < iend - istart; ++i) { - // THIS LOOKS LIKE IT COULD INTRODUCE A SHARP DISCONTINUITY. COULD THAT BE A SOURCE OF NUMERICAL PROBLEMS? MAYBE, BUT NOTHING TOO SIGNIFICANT THOUGH - const PetscScalar norm = - BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); - snes_x_data[i] /= norm; // FIXME: Is this right? In most cases this means the variable is just 1. Yes, because this is just the guess for the next solution - x1_data[i] /= norm; // Update history for predictor - var_scaling_factors_data[i] *= norm; - } + // Individual variable scaling + // Note: If variables are rescaled then the Jacobian columns + // need to be scaled or recalculated + int istart, iend; + VecGetOwnershipRange(snes_x, &istart, &iend); + + // Take ownership of snes_x and var_scaling_factors data + PetscScalar* snes_x_data = nullptr; + PetscCall(VecGetArray(snes_x, &snes_x_data)); + PetscScalar* x1_data; + PetscCall(VecGetArray(x1, &x1_data)); + PetscScalar* var_scaling_factors_data; + PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); + + // Normalise each value in the state + // Limit normalisation so scaling factor is never smaller than rtol + for (int i = 0; i < iend - istart; ++i) { + const PetscScalar norm = + BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); + snes_x_data[i] /= + norm; + x1_data[i] /= norm; // Update history for predictor + var_scaling_factors_data[i] *= norm; + } - // Restore vector underlying data - PetscCall(VecRestoreArray(var_scaling_factors, &var_scaling_factors_data)); - PetscCall(VecRestoreArray(x1, &x1_data)); - PetscCall(VecRestoreArray(snes_x, &snes_x_data)); + // Restore vector underlying data + PetscCall(VecRestoreArray(var_scaling_factors, &var_scaling_factors_data)); + PetscCall(VecRestoreArray(x1, &x1_data)); + PetscCall(VecRestoreArray(snes_x, &snes_x_data)); - if (diagnose) { - // Print maximum and minimum scaling factors - PetscReal max_scale, min_scale; - VecMax(var_scaling_factors, nullptr, &max_scale); - VecMin(var_scaling_factors, nullptr, &min_scale); - output.write("Var scaling: {} -> {}\n", min_scale, max_scale); - } - - // Force recalculation of the Jacobian - // TODO: Get my head around lagged jacobians so I can check if this is right - // SNESGetLagJacobianPersists(snes, saved_jacobian_persist); - // SNESSetLagJacobianPersists(snes, PETSC_FALSE); + if (diagnose) { + // Print maximum and minimum scaling factors + PetscReal max_scale, min_scale; + VecMax(var_scaling_factors, nullptr, &max_scale); + VecMin(var_scaling_factors, nullptr, &min_scale); + output.write("Var scaling: {} -> {}\n", min_scale, max_scale); + } - SNESGetLagJacobian(snes, &saved_jacobian_lag); - SNESSetLagJacobian(snes, 1); - return PETSC_SUCCESS; + // Force recalculation of the Jacobian + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + return PETSC_SUCCESS; } int SNESSolver::run() { @@ -965,12 +962,12 @@ int SNESSolver::run() { // If saving solver residuals, need to get initial value. if (diagnose) { snes_function(snes_x, snes_f, false); - } + } BoutReal target = simtime; - BoutReal accumulated_change = 0.; recent_failure_rate = 0.0; - int loop_count = 0; + int steps_since_rescale = 0; + BoutReal change_since_rescale = 0.; for (int s = 1; s <= getNumberOutputSteps(); s++) { target += getOutputTimestep(); @@ -987,11 +984,12 @@ int SNESSolver::run() { break; // Could happen if step over multiple outputs } - //output.write("Accumulated change to state is: {}", accumulated_change); - if (scale_vars and (accumulated_change > rescale_threshold or loop_count == rescale_period)) { + if (scale_vars + and (change_since_rescale > rescale_threshold + or steps_since_rescale == rescale_period)) { PetscCall(rescale(saved_jacobian_lag)); - accumulated_change = 0.; - loop_count = 0; + change_since_rescale = 0.; + steps_since_rescale = 0; } // Copy the state (snes_x) into initial values (x0) @@ -1037,8 +1035,6 @@ int SNESSolver::run() { // Note: When the timestep is changed the preconditioner needs to be updated // => Step over the output time and interpolate if not matrix free - // FIXME: Is this the most appropriate thing to do when we're taking such massive timesteps? - if (matrix_free) { // Ensure that the timestep goes to the next output time and then stops. // This avoids the need to interpolate @@ -1126,9 +1122,9 @@ int SNESSolver::run() { // Restore state VecCopy(x0, snes_x); // FIXME: Not sure I need this - if (diagnose) { - VecCopy(f0, snes_f); - } + // if (diagnose) { + // VecCopy(f0, snes_f); + // } // Recalculate the Jacobian if (jacobian_pruned and (snes_failures > 2) and (4 * lin_its > 3 * maxl)) { @@ -1140,7 +1136,6 @@ int SNESSolver::run() { PetscCall(FDJrestoreFromPruning()); } - // **** FIXME: Replace with booleans if (saved_jacobian_lag == 0) { // This triggers a Jacobian recalculation SNESGetLagJacobian(snes, &saved_jacobian_lag); @@ -1161,7 +1156,6 @@ int SNESSolver::run() { continue; // Try again } - // **** FIXME: Replace with booleans if (saved_jacobian_lag != 0) { // Following successful step, reset Jacobian lag // to previous value @@ -1213,8 +1207,7 @@ int SNESSolver::run() { PetscCall(VecRestoreArray(deriv, &fdata)); // Forward Euler VecAXPY(snes_x, dt, deriv); - } - else { + } else { BoutReal* fdata = nullptr; PetscCall(VecGetArray(snes_f, &fdata)); save_derivs(fdata); @@ -1222,24 +1215,21 @@ int SNESSolver::run() { // Forward Euler VecAXPY(snes_x, dt, snes_f); } - } - // int i; - // PetscReal norm, minval, maxval; - // PetscCall(VecNorm(snes_f, NORM_2, &norm)); - // PetscCall(VecMin(snes_f, &i, &minval)); - // PetscCall(VecMax(snes_f, &i, &maxval)); - // output.write("Max, min, and norm of residual for SNES function: {}, {}, {}\n", maxval, minval, norm); - + // int i; + // PetscReal norm, minval, maxval; + // PetscCall(VecNorm(snes_f, NORM_2, &norm)); + // PetscCall(VecMin(snes_f, &i, &minval)); + // PetscCall(VecMax(snes_f, &i, &maxval)); + // output.write("Max, min, and norm of residual for SNES function: {}, {}, {}\n", maxval, minval, norm); simtime += dt; // Update local and global residuals - /// FIXME: This overwrites snes_f with the time derivs! PetscCall(updateResiduals(snes_x)); - accumulated_change += global_residual * dt; - ++loop_count; + change_since_rescale += global_residual * dt; + ++steps_since_rescale; if (diagnose) { // Gather and print diagnostic information @@ -1307,7 +1297,6 @@ int SNESSolver::run() { VecAXPBY(output_x, alpha, 1. - alpha, x0); if (diagnose) { - // FIXME: Is this a good approach to take here? It's not the *actual* residual of a solve this way. VecCopy(snes_f, output_f); VecAXPBY(output_f, alpha, 1 - alpha, f0); } @@ -1326,7 +1315,6 @@ int SNESSolver::run() { // scaled_x <- output_x * var_scaling_factors PetscCall(VecPointwiseMult(scaled_x, output_x, var_scaling_factors)); } else if (asinh_vars) { - // FIXME: scale_vars and asinh_vars aren't mutually exclusive, are they? PetscCall(VecCopy(output_x, scaled_x)); } else { scaled_x = output_x; @@ -1422,17 +1410,17 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { const BoutReal* current_residual = nullptr; if (diagnose) { - // Call RHS function to get time derivatives + // Call RHS function to get time derivatives PetscCall(rhs_function(x, deriv, false)); - // Reading the residual vectors - PetscCall(VecGetArrayRead(deriv, ¤t_residual)); + // Reading the residual vectors + PetscCall(VecGetArrayRead(deriv, ¤t_residual)); } else { - // Call RHS function to get time derivatives - PetscCall(rhs_function(x, snes_f, false)); + // Call RHS function to get time derivatives + PetscCall(rhs_function(x, snes_f, false)); - // Reading the residual vectors - PetscCall(VecGetArrayRead(snes_f, ¤t_residual)); + // Reading the residual vectors + PetscCall(VecGetArrayRead(snes_f, ¤t_residual)); } // Note: The ordering of quantities in the PETSc vectors @@ -1502,7 +1490,7 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { } if (diagnose) { - // Restore Vec data arrays + // Restore Vec data arrays PetscCall(VecRestoreArrayRead(deriv, ¤t_residual)); } else { PetscCall(VecRestoreArrayRead(snes_f, ¤t_residual)); diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 82f070805a..9b5a6a6eeb 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -206,8 +206,8 @@ private: BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes) 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 - BoutReal last_failure_weight; ///< 1 / number of recent solves + BoutReal recent_failure_rate; ///< Rolling average of recent failure rate + BoutReal last_failure_weight; ///< 1 / number of recent solves BoutReal nl_its_prev; BoutReal nl_its_prev2; @@ -222,10 +222,10 @@ private: PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function - Vec deriv; ///< Time derivative; only used if diagnose = true, otherwise will store in snes_f + Vec deriv; ///< Time derivative; only used if diagnose = true, otherwise will store in snes_f Vec snes_x; ///< Result of SNES Vec x0; ///< Solution at start of current timestep - Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) + Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) Vec delta_x; ///< Change in solution Vec output_x; ///< Solution to output. Used if interpolating. Vec output_f; ///< Residual to output, if diagnose == true. Used if interpolating. @@ -268,15 +268,18 @@ private: bool scale_vars; ///< Scale individual variables? int rescale_period; ///< How many time-steps before rescaling variables - BoutReal rescale_threshold; //< How much change in the state there should be before rescaling + BoutReal + rescale_threshold; //< How much change in the state there should be before rescaling Vec var_scaling_factors; ///< Factors to multiply variables when passing to user 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 - std::vector resid_2d; ///< Storage for residuals of SNES solve, unpacked from snes_f - std::vector resid_3d; ///< Storage for residuals of SNES solve, unpacked from snes_f + std::vector + resid_2d; ///< Storage for residuals of SNES solve, unpacked from snes_f + std::vector + resid_3d; ///< Storage for residuals of SNES solve, unpacked from snes_f }; #else diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 0afb21f528..b46b890574 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1038,7 +1038,8 @@ void Solver::load_vars(BoutReal* udata) { f.var->setLocation(f.location); } - loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); + loop_vars(VarRange(f2d), + VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); // Mark each vector as either co- or contra-variant @@ -1060,7 +1061,8 @@ void Solver::load_derivs(BoutReal* udata) { f.F_var->setLocation(f.location); } - loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); + loop_vars(VarRange(f2d), + VarRange(f3d), udata, SOLVER_VAR_OP::LOAD); // Mark each vector as either co- or contra-variant @@ -1102,7 +1104,8 @@ void Solver::save_vars(BoutReal* udata) { } } - loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::SAVE); + loop_vars(VarRange(f2d), + VarRange(f3d), udata, SOLVER_VAR_OP::SAVE); } void Solver::save_derivs(BoutReal* dudata) { @@ -1132,10 +1135,14 @@ void Solver::save_derivs(BoutReal* dudata) { } } - loop_vars(VarRange(f2d), VarRange(f3d), dudata, SOLVER_VAR_OP::SAVE); + loop_vars(VarRange(f2d), + VarRange(f3d), dudata, SOLVER_VAR_OP::SAVE); } -void Solver::set_id(BoutReal* udata) { loop_vars(VarRange(f2d), VarRange(f3d), udata, SOLVER_VAR_OP::SET_ID); } +void Solver::set_id(BoutReal* udata) { + loop_vars(VarRange(f2d), + VarRange(f3d), udata, SOLVER_VAR_OP::SET_ID); +} Field3D Solver::globalIndex(int localStart) { // Use global mesh: FIX THIS! From 87d29ff2095874805b75d7c14375ecc25db3169c Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 19 Jun 2026 17:56:11 +0100 Subject: [PATCH 06/13] Add documentation for new approach to scaling in SNES --- manual/sphinx/user_docs/time_integration.rst | 32 ++++++++++++++++++++ src/solver/impls/snes/snes.cxx | 2 +- 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index a2bd896f85..aae0d4fe9d 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -747,6 +747,34 @@ Setting ``solver:force_symmetric_coloring = true``, will make sure that the jacobian colouring matrix is symmetric. This will often include a few extra non-zeros that the stencil will miss otherwise + +Variable Scaling +~~~~~~~~~~~~~~~~ + +There may be differences of many orders of magnitude between your +variables or within variables across the domain. This can result in a +particular area of the domain for a particular variable dominating the +residual in the nonlinear solve because its residual has the largest +absolute value, even if not the largest relative value. As a +consequence, tighter tolerances will be needed to ensure other +variables and parts of the domain are solved to sufficient +accuracy. The ``scale_vars`` option can help address this by +renormalising all variables to be of order unity across the entire +domain. + +.. code-block:: ini + + scale_vars = true + rescale_period = 30 # Maximum number of time-steps taken before rescaling the variables + rescale_threshold = 100. # Approximate overall change to variables permitted before rescaling + +It has been found that scaling variables in this way allows +simulations to run with much looser tolerances than would otherwise be +possible (e.g., ``rtol = 1e-5`` and ``atol = 1e-3``). Four-times +speedups have been observed by doing this. Once a steady-state has +been reached the simulation can be run for a further few time-steps +with tighter tolerances to improve the accuracy. + Diagnostics and Monitoring --------------------------- @@ -765,6 +793,10 @@ When ``equation_form = pseudo_transient``, the solver saves additional diagnosti These can be visualized to understand convergence behavior and identify problematic regions. +The residuals from the last nonlinear solve are also saved with names +``resid_``. Plotting these can help to understand which +variables and parts of the domain are controlling convergence. + Summary of solver options ~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 21757cda90..571dc9504e 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -630,7 +630,7 @@ SNESSolver::SNESSolver(Options* opts) .doc("Scale time derivatives (Jacobian row scaling)?") .withDefault(false)), scale_vars((*options)["scale_vars"] - .doc("Scale variables (Jacobian column scaling)?") + .doc("Scale variables to be order unity?") .withDefault(false)), rescale_period( (*options)["rescale_period"] From 81fda045973e9e9877cbad8016ef52f68c1b3ac0 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 19 Jun 2026 20:02:34 +0100 Subject: [PATCH 07/13] Fixed linting complaints --- include/bout/solver.hxx | 13 ++++++++----- src/solver/impls/snes/snes.cxx | 25 ++++++++++++++----------- src/solver/impls/snes/snes.hxx | 16 +++++++++------- src/solver/solver.cxx | 2 +- 4 files changed, 32 insertions(+), 24 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index fa67ac1ca7..eba8164419 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -40,10 +40,13 @@ #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" +#include "bout/globals.hxx" #include "bout/monitor.hxx" #include "bout/options.hxx" +#include "bout/regions.hxx" #include "bout/unused.hxx" +#include #include /////////////////////////////////////////////////////////////////// @@ -92,8 +95,8 @@ constexpr auto SOLVERIMEXBDF2 = "imexbdf2"; constexpr auto SOLVERSNES = "snes"; constexpr auto SOLVERRKGENERIC = "rkgeneric"; -enum class FieldCategories { VARS, DERIVS, MMS }; -enum class SOLVER_VAR_OP { LOAD, SET_ID, SAVE }; +enum class FieldCategories : std::uint8_t { VARS, DERIVS, MMS }; +enum class SOLVER_VAR_OP : std::uint8_t { LOAD, SET_ID, SAVE }; /// A type to set where in the list monitors are added enum class MonitorPosition { BACK, FRONT }; @@ -635,7 +638,7 @@ protected: void loop_vars(RangeF2D f2d_range, RangeF3D f3d_range, BoutReal* udata, SOLVER_VAR_OP op) { // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; + const Mesh* mesh = bout::globals::mesh; int p = 0; // Counter for location in udata array @@ -719,9 +722,9 @@ private: * Is the interleaving of variables needed or helpful to the solver? **************************************************************************/ // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; + const Mesh* mesh = bout::globals::mesh; - int nz = mesh->LocalNz; + const int nz = mesh->LocalNz; switch (op) { case SOLVER_VAR_OP::LOAD: { diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 571dc9504e..3c26db4d5b 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -885,15 +885,16 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { // Individual variable scaling // Note: If variables are rescaled then the Jacobian columns // need to be scaled or recalculated - int istart, iend; + int istart = 0; + int iend = 0; VecGetOwnershipRange(snes_x, &istart, &iend); // Take ownership of snes_x and var_scaling_factors data PetscScalar* snes_x_data = nullptr; PetscCall(VecGetArray(snes_x, &snes_x_data)); - PetscScalar* x1_data; + PetscScalar* x1_data = nullptr; PetscCall(VecGetArray(x1, &x1_data)); - PetscScalar* var_scaling_factors_data; + PetscScalar* var_scaling_factors_data = nullptr; PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); // Normalise each value in the state @@ -901,8 +902,7 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { for (int i = 0; i < iend - istart; ++i) { const PetscScalar norm = BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); - snes_x_data[i] /= - norm; + snes_x_data[i] /= norm; x1_data[i] /= norm; // Update history for predictor var_scaling_factors_data[i] *= norm; } @@ -914,7 +914,8 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { if (diagnose) { // Print maximum and minimum scaling factors - PetscReal max_scale, min_scale; + PetscReal max_scale = 0.; + PetscReal min_scale = 0.; VecMax(var_scaling_factors, nullptr, &max_scale); VecMin(var_scaling_factors, nullptr, &min_scale); output.write("Var scaling: {} -> {}\n", min_scale, max_scale); @@ -922,6 +923,10 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { // Force recalculation of the Jacobian SNESGetLagJacobian(snes, &saved_jacobian_lag); + // FIXME: This isn't actually what we want to happen. The lag should + // stay as before, we just want to make sure the Jacobian gets + // re-evaluated at the start of the solve. We need to use + // SNESSetJacobianPersist for that. SNESSetLagJacobian(snes, 1); return PETSC_SUCCESS; } @@ -945,8 +950,8 @@ int SNESSolver::run() { PetscCall(VecRestoreArray(snes_x, &xdata)); } + int saved_jacobian_lag = 0; if (scale_vars) { - int saved_jacobian_lag; PetscCall(rescale(saved_jacobian_lag)); } @@ -973,8 +978,6 @@ int SNESSolver::run() { bool looping = true; int snes_failures = 0; // Count SNES convergence failures - int saved_jacobian_lag = 0; - PetscBool saved_jacobian_persist = PETSC_FALSE; const BoutReal start_global_residual = global_residual; do { @@ -1347,7 +1350,7 @@ int SNESSolver::run() { } if (diagnose) { - BoutReal* resid_data; + BoutReal* resid_data = nullptr; PetscCall(VecGetArray(output_f, &resid_data)); loop_vars(resid_2d, resid_3d, resid_data, SOLVER_VAR_OP::LOAD); PetscCall(VecRestoreArray(output_f, &resid_data)); @@ -1789,7 +1792,7 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // if (scale_vars) { // VecPointwiseMult(delta_x, delta_x, var_scaling_factors); // } - VecAXPY(f, -1. / dt, delta_x); // f <- f - delta_x / dt + VecAXPY(f, -1. / dt, delta_x); // f <- f - delta_x / dt break; } case BoutSnesEquationForm::pseudo_transient: { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 9b5a6a6eeb..3b1ffe99f0 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -35,6 +35,8 @@ class SNESSolver; +#include + #include "mpi.h" #include @@ -223,10 +225,10 @@ private: PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function Vec deriv; ///< Time derivative; only used if diagnose = true, otherwise will store in snes_f - Vec snes_x; ///< Result of SNES - Vec x0; ///< Solution at start of current timestep - Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) - Vec delta_x; ///< Change in solution + Vec snes_x; ///< Result of SNES + Vec x0; ///< Solution at start of current timestep + Vec f0; ///< Residual at start of current timestep (only stored if diagnose = true) + Vec delta_x; ///< Change in solution Vec output_x; ///< Solution to output. Used if interpolating. Vec output_f; ///< Residual to output, if diagnose == true. Used if interpolating. @@ -252,7 +254,7 @@ private: bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? int lag_jacobian; ///< Re-use Jacobian bool jacobian_persists; ///< Re-use Jacobian and preconditioner across nonlinear solves - bool use_coloring; ///< Use matrix coloring + bool use_coloring; ///< Use matrix coloring bool jacobian_recalculated; ///< Flag set when Jacobian is recalculated bool prune_jacobian; ///< Remove small elements in the Jacobian? @@ -266,8 +268,8 @@ private: Vec rhs_scaling_factors; ///< Factors to multiply RHS function Vec jac_row_inv_norms; ///< 1 / Norm of the rows of the Jacobian - bool scale_vars; ///< Scale individual variables? - int rescale_period; ///< How many time-steps before rescaling variables + bool scale_vars; ///< Scale individual variables? + int rescale_period; ///< How many time-steps before rescaling variables BoutReal rescale_threshold; //< How much change in the state there should be before rescaling Vec var_scaling_factors; ///< Factors to multiply variables when passing to user diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index b46b890574..1f7b654427 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -2,7 +2,7 @@ * Copyright 2010 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify From 70d54848cc92d04cfb2a18d989522602836edc28 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 19 Jun 2026 20:46:58 +0100 Subject: [PATCH 08/13] Fixed a few bugs and clang-tidy issues --- include/bout/solver.hxx | 6 ++++-- src/solver/impls/snes/snes.cxx | 8 ++++---- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index eba8164419..1a3d475edf 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -41,11 +41,13 @@ #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/globals.hxx" +#include "bout/mesh.hxx" #include "bout/monitor.hxx" #include "bout/options.hxx" -#include "bout/regions.hxx" +#include "bout/region.hxx" #include "bout/unused.hxx" +#include #include #include @@ -487,7 +489,7 @@ protected: } private: - underlying_iterator it; + underlying_iterator it{}; }; template diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 3c26db4d5b..5003117e3f 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1204,10 +1204,10 @@ int SNESSolver::run() { // Copy derivatives back if (diagnose) { - BoutReal* fdata = nullptr; - PetscCall(VecGetArray(deriv, &fdata)); - save_derivs(fdata); - PetscCall(VecRestoreArray(deriv, &fdata)); + BoutReal* ddata = nullptr; + PetscCall(VecGetArray(deriv, &ddata)); + save_derivs(ddata); + PetscCall(VecRestoreArray(deriv, &ddata)); // Forward Euler VecAXPY(snes_x, dt, deriv); } else { From 1a66918f5caf48b838e78d058d6c694378792273 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 20 Jun 2026 13:14:42 -0700 Subject: [PATCH 09/13] Solver::VarIterator: Replace switch with if constexpr C is a template parameter so can be determined at compile time. --- include/bout/solver.hxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 1a3d475edf..24398c8fd2 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -410,22 +410,22 @@ protected: VarIterator(underlying_iterator _it) : it(std::move(_it)) {} reference operator*() const { - switch (C) { - case (FieldCategories::VARS): + if constexpr (C == FieldCategories::VARS) { return *(it->var); - case (FieldCategories::DERIVS): + } else if constexpr (C == FieldCategories::DERIVS) { return *(it->F_var); - case (FieldCategories::MMS): + } else { + static_assert(C == FieldCategories::MMS); return *(it->MMS_err); } } pointer operator->() const { - switch (C) { - case (FieldCategories::VARS): + if constexpr (C == FieldCategories::VARS) { return it->var; - case (FieldCategories::DERIVS): + } else if constexpr (C == FieldCategories::DERIVS) { return it->F_var; - case (FieldCategories::MMS): + } else { + static_assert(C == FieldCategories::MMS); return it->MMS_err.get(); } } From af78c6f293100b7f71cb8916eb6688826662e4d2 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 20 Jun 2026 13:19:56 -0700 Subject: [PATCH 10/13] SNESSolver::rescale: Handle predictor=false case The x1 vector is only allocated if predictor is enabled. --- src/solver/impls/snes/snes.cxx | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 5003117e3f..db34eb43f7 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -893,7 +893,10 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { PetscScalar* snes_x_data = nullptr; PetscCall(VecGetArray(snes_x, &snes_x_data)); PetscScalar* x1_data = nullptr; - PetscCall(VecGetArray(x1, &x1_data)); + if (predictor) { + // x1 is only allocated if predictor is enabled + PetscCall(VecGetArray(x1, &x1_data)); + } PetscScalar* var_scaling_factors_data = nullptr; PetscCall(VecGetArray(var_scaling_factors, &var_scaling_factors_data)); @@ -903,13 +906,17 @@ PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) { const PetscScalar norm = BOUTMAX(std::abs(snes_x_data[i]), rtol / var_scaling_factors_data[i]); snes_x_data[i] /= norm; - x1_data[i] /= norm; // Update history for predictor + if (predictor) { + x1_data[i] /= norm; // Update history for predictor + } var_scaling_factors_data[i] *= norm; } // Restore vector underlying data PetscCall(VecRestoreArray(var_scaling_factors, &var_scaling_factors_data)); - PetscCall(VecRestoreArray(x1, &x1_data)); + if (predictor) { + PetscCall(VecRestoreArray(x1, &x1_data)); + } PetscCall(VecRestoreArray(snes_x, &snes_x_data)); if (diagnose) { From 2cfb90a582646d9acb9a275ef332a4645fc6cbb6 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 20 Jun 2026 13:26:32 -0700 Subject: [PATCH 11/13] Solver::loop_vars_op: Increment metadata unconditionally Previously metadata_it only advanced on the non-continue path. If one field had evolve_bndry == false, every later field in that boundary loop would be checked against the same stale metadata entry and skipped too. --- include/bout/solver.hxx | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 24398c8fd2..7030559cd9 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -735,12 +735,13 @@ private: // Loop over 2D variables auto metadata_it = f2d.begin(); for (auto& field : f2d_range) { - if (bndry && !metadata_it->evolve_bndry) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { continue; } field[i2d] = udata[p]; p++; - metadata_it++; } for (int jz = 0; jz < nz; jz++) { @@ -748,12 +749,13 @@ private: auto metadata_it = f3d.begin(); // Loop over 3D variables for (auto& field : f3d_range) { - if (bndry && !metadata_it->evolve_bndry) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { continue; } field[field.getMesh()->ind2Dto3D(i2d, jz)] = udata[p]; p++; - metadata_it++; } } break; @@ -798,12 +800,13 @@ private: // Loop over 2D variables auto metadata_it = f2d.begin(); for (const auto& field : f2d_range) { - if (bndry && !metadata_it->evolve_bndry) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { continue; } udata[p] = field[i2d]; p++; - metadata_it++; } for (int jz = 0; jz < nz; jz++) { @@ -811,12 +814,13 @@ private: auto metadata_it = f3d.begin(); // Loop over 3D variables for (const auto& field : f3d_range) { - if (bndry && !metadata_it->evolve_bndry) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { continue; } udata[p] = field[field.getMesh()->ind2Dto3D(i2d, jz)]; p++; - metadata_it++; } } break; From 013aef48f8e1f464686f265ae9c9c2a701b25345 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 20 Jun 2026 13:30:07 -0700 Subject: [PATCH 12/13] SNESSolver: scale_vars in zero iteration case When SNES takes zero iterations (nl_its == 0) a simple Euler step is taken. If scale_vars is true then the derivatives should be scaled to be consistent with rhs_function. --- src/solver/impls/snes/snes.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index db34eb43f7..5f3896a5c9 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1215,6 +1215,9 @@ int SNESSolver::run() { PetscCall(VecGetArray(deriv, &ddata)); save_derivs(ddata); PetscCall(VecRestoreArray(deriv, &ddata)); + if (scale_vars) { + PetscCall(VecPointwiseDivide(deriv, deriv, var_scaling_factors)); + } // Forward Euler VecAXPY(snes_x, dt, deriv); } else { @@ -1222,6 +1225,9 @@ int SNESSolver::run() { PetscCall(VecGetArray(snes_f, &fdata)); save_derivs(fdata); PetscCall(VecRestoreArray(snes_f, &fdata)); + if (scale_vars) { + PetscCall(VecPointwiseDivide(snes_f, snes_f, var_scaling_factors)); + } // Forward Euler VecAXPY(snes_x, dt, snes_f); } From e0ca819e2fd5d0b7c17110826c011392b6d095a2 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 20 Jun 2026 14:02:43 -0700 Subject: [PATCH 13/13] Delete commented-out code --- include/bout/solver.hxx | 10 +++++----- src/solver/impls/snes/snes.cxx | 23 +---------------------- 2 files changed, 6 insertions(+), 27 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 7030559cd9..09ede6a32b 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -718,11 +718,11 @@ private: void loop_vars_op(RangeF2D f2d_range, RangeF3D f3d_range, Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, bool bndry) { /************************************************************************** - * Looping over variables - * - * NOTE: This part is very inefficient, and should be replaced ASAP - * Is the interleaving of variables needed or helpful to the solver? - **************************************************************************/ + * Looping over variables + * + * NOTE: This part is very inefficient, and should be replaced ASAP + * Is the interleaving of variables needed or helpful to the solver? + **************************************************************************/ // Use global mesh: FIX THIS! const Mesh* mesh = bout::globals::mesh; diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 5f3896a5c9..b703ab1c90 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1131,10 +1131,6 @@ int SNESSolver::run() { } // Restore state VecCopy(x0, snes_x); - // FIXME: Not sure I need this - // if (diagnose) { - // VecCopy(f0, snes_f); - // } // Recalculate the Jacobian if (jacobian_pruned and (snes_failures > 2) and (4 * lin_its > 3 * maxl)) { @@ -1233,13 +1229,6 @@ int SNESSolver::run() { } } - // int i; - // PetscReal norm, minval, maxval; - // PetscCall(VecNorm(snes_f, NORM_2, &norm)); - // PetscCall(VecMin(snes_f, &i, &minval)); - // PetscCall(VecMax(snes_f, &i, &maxval)); - // output.write("Max, min, and norm of residual for SNES function: {}, {}, {}\n", maxval, minval, norm); - simtime += dt; // Update local and global residuals @@ -1802,10 +1791,7 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // f = (x0 - x)/Δt + f // First calculate x - x0 to minimise floating point issues VecWAXPY(delta_x, -1.0, x0, x); // delta_x = x - x0 - // if (scale_vars) { - // VecPointwiseMult(delta_x, delta_x, var_scaling_factors); - // } - VecAXPY(f, -1. / dt, delta_x); // f <- f - delta_x / dt + VecAXPY(f, -1. / dt, delta_x); // f <- f - delta_x / dt break; } case BoutSnesEquationForm::pseudo_transient: { @@ -1813,9 +1799,6 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // except that Δt is a vector // f = (x0 - x)/Δt + f VecWAXPY(delta_x, -1.0, x0, x); - // if (scale_vars) { - // VecPointwiseMult(delta_x, delta_x, var_scaling_factors); - // } VecPointwiseDivide(delta_x, delta_x, dt_vec); // delta_x /= dt VecAXPY(f, -1., delta_x); // f <- f - delta_x break; @@ -1823,10 +1806,6 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { case BoutSnesEquationForm::backward_euler: { // Backward Euler // Set f = x - x0 - Δt*f - // if (scale_vars) { - // VecPointwiseMult(x, x, var_scaling_factors); - // VecPointwiseMult(x0, x0, var_scaling_factors); - // } VecAYPX(f, -dt, x); // f <- x - Δt*f VecAXPY(f, -1.0, x0); // f <- f - x0 break;