diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 2e0da05224..09ede6a32b 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -40,10 +40,15 @@ #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/region.hxx" #include "bout/unused.hxx" +#include +#include #include /////////////////////////////////////////////////////////////////// @@ -92,7 +97,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 : 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 }; @@ -391,6 +397,114 @@ 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 { + if constexpr (C == FieldCategories::VARS) { + return *(it->var); + } else if constexpr (C == FieldCategories::DERIVS) { + return *(it->F_var); + } else { + static_assert(C == FieldCategories::MMS); + return *(it->MMS_err); + } + } + pointer operator->() const { + if constexpr (C == FieldCategories::VARS) { + return it->var; + } else if constexpr (C == FieldCategories::DERIVS) { + return it->F_var; + } else { + static_assert(C == 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 +633,28 @@ 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! + const 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 +709,124 @@ 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! + const Mesh* mesh = bout::globals::mesh; + + const 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) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { + continue; + } + field[i2d] = udata[p]; + p++; + } + + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (auto& field : f3d_range) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { + continue; + } + field[field.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& 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) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { + continue; + } + udata[p] = field[i2d]; + p++; + } + + for (int jz = 0; jz < nz; jz++) { + + auto metadata_it = f3d.begin(); + // Loop over 3D variables + for (const auto& field : f3d_range) { + const auto& metadata = *metadata_it; + ++metadata_it; + if (bndry && !metadata.evolve_bndry) { + continue; + } + udata[p] = field[field.getMesh()->ind2Dto3D(i2d, jz)]; + p++; + } + } + break; + } + } + } /// Check if a variable has already been added bool varAdded(const std::string& name); 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 1dd5e37c84..b703ab1c90 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -630,8 +630,16 @@ 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"] + .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)) {} @@ -655,6 +663,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 +683,10 @@ int SNESSolver::init() { PetscCall(VecDuplicate(snes_x, &snes_f)); PetscCall(VecDuplicate(snes_x, &x0)); + if (diagnose) { + PetscCall(VecDuplicate(snes_f, &f0)); + PetscCall(VecDuplicate(snes_f, &deriv)); + } if ((equation_form == BoutSnesEquationForm::rearranged_backward_euler) || (equation_form == BoutSnesEquationForm::pseudo_transient)) { @@ -755,6 +777,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 +881,63 @@ 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 = 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 = nullptr; + 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)); + + // 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; + 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)); + if (predictor) { + PetscCall(VecRestoreArray(x1, &x1_data)); + } + PetscCall(VecRestoreArray(snes_x, &snes_x_data)); + + if (diagnose) { + // Print maximum and minimum scaling factors + 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); + } + + // 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; +} + int SNESSolver::run() { // Set initial guess at the solution from variables { @@ -875,6 +957,11 @@ int SNESSolver::run() { PetscCall(VecRestoreArray(snes_x, &xdata)); } + int saved_jacobian_lag = 0; + if (scale_vars) { + PetscCall(rescale(saved_jacobian_lag)); + } + // Initialise residuals local_residual = 0.0; local_residual_2d = 0.0; @@ -884,15 +971,20 @@ 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; + int steps_since_rescale = 0; + BoutReal change_since_rescale = 0.; for (int s = 1; s <= getNumberOutputSteps(); s++) { target += getOutputTimestep(); bool looping = true; int snes_failures = 0; // Count SNES convergence failures - int saved_jacobian_lag = 0; - int loop_count = 0; const BoutReal start_global_residual = global_residual; do { @@ -902,57 +994,19 @@ 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 (change_since_rescale > rescale_threshold + or steps_since_rescale == rescale_period)) { + PetscCall(rescale(saved_jacobian_lag)); + change_since_rescale = 0.; + steps_since_rescale = 0; } - ++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 @@ -1152,21 +1206,35 @@ int SNESSolver::run() { } // Copy derivatives back - { + if (diagnose) { + BoutReal* ddata = nullptr; + 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 { BoutReal* fdata = nullptr; 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); } - - // Forward Euler - VecAXPY(snes_x, dt, snes_f); } simtime += dt; // Update local and global residuals PetscCall(updateResiduals(snes_x)); + change_since_rescale += global_residual * dt; + ++steps_since_rescale; if (diagnose) { // Gather and print diagnostic information @@ -1232,11 +1300,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 @@ -1275,6 +1351,13 @@ int SNESSolver::run() { BoutComm::abort(1); } + if (diagnose) { + 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)); + } + if (call_monitors(output_time, s, getNumberOutputSteps()) != 0) { break; // User signalled to quit } @@ -1330,12 +1413,20 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { local_residual_2d_prev = copy(local_residual_2d); global_residual_prev = global_residual; - // 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)); + 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 + PetscCall(VecGetArrayRead(snes_f, ¤t_residual)); + } // Note: The ordering of quantities in the PETSc vectors // depends on the Solver::loop_vars function @@ -1403,8 +1494,12 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { } } - // Restore Vec data arrays - PetscCall(VecRestoreArrayRead(snes_f, ¤t_residual)); + if (diagnose) { + // Restore Vec data arrays + 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)); @@ -1672,6 +1767,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; } @@ -1933,6 +2032,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..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 @@ -115,6 +117,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 @@ -202,8 +208,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; @@ -218,10 +224,13 @@ private: PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function - Vec snes_x; ///< Result of SNES - Vec x0; ///< Solution at start of current timestep - Vec delta_x; ///< Change in solution + 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 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 @@ -245,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? @@ -259,12 +268,20 @@ 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? + 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 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..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 @@ -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,8 @@ 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 +1061,8 @@ 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 +1104,8 @@ 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 +1135,14 @@ 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!