From 0ced35ef3c911b77a77e9d1ca0247bdc28ebc30e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 2 May 2025 21:10:28 -0700 Subject: [PATCH 1/2] beuler: Add output interpolation The preconditioner depends on the timestep, so avoid changing timestep unless necessary. Allow timesteps to go past output time, then interpolate. If timesteps are adjusted then recalculate the Jacobian. This could be optimised to update the Jacobian and recalculate the preconditioner, without finite differencing the RHS. --- src/solver/impls/snes/snes.cxx | 96 ++++++++++++++++++++++++++++------ src/solver/impls/snes/snes.hxx | 17 +++--- 2 files changed, 89 insertions(+), 24 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 4dce1e7086..0e514d56cc 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -297,13 +297,21 @@ int SNESSolver::init() { SNESSetJacobian(snes, Jmf, Jmf, MatMFFDComputeJacobian, this); } else { - // Calculate the Jacobian using finite differences. - // The finite difference Jacobian (Jfd) may be used for both operator - // and preconditioner or, if matrix_free_operator, in only the preconditioner. + // Calculate the Jacobian using finite differences. The finite + // difference Jacobian (Jfd) may be used for both operator and + // preconditioner or, if matrix_free_operator, in only the + // preconditioner. + + // Create a vector to store interpolated output solution + // Used so that the timestep does not have to be adjusted, + // because that would require updating the preconditioner. + ierr = VecDuplicate(snes_x, &output_x); + CHKERRQ(ierr); + if (use_coloring) { - // Use matrix coloring - // This greatly reduces the number of times the rhs() function needs - // to be evaluated when calculating the Jacobian. + // Use matrix coloring. + // This greatly reduces the number of times the rhs() function + // needs to be evaluated when calculating the Jacobian. // Use global mesh for now Mesh* mesh = bout::globals::mesh; @@ -717,14 +725,17 @@ int SNESSolver::run() { CHKERRQ(ierr); } + BoutReal target = simtime; for (int s = 0; s < getNumberOutputSteps(); s++) { - BoutReal target = simtime + getOutputTimestep(); + target += getOutputTimestep(); bool looping = true; int snes_failures = 0; // Count SNES convergence failures int saved_jacobian_lag = 0; int loop_count = 0; do { + if (simtime >= target) + break; // Could happen if step over multiple outputs if (scale_vars) { // Individual variable scaling // Note: If variables are rescaled then the Jacobian columns @@ -802,9 +813,15 @@ int SNESSolver::run() { dt = timestep; looping = true; if (simtime + dt >= target) { - // Ensure that the timestep goes to the next output time and then stops + // Note: When the timestep is changed the preconditioner needs to be updated + // => Step over the output time and interpolate if not matrix free + + if (matrix_free) { + // Ensure that the timestep goes to the next output time and then stops. + // This avoids the need to interpolate + dt = target - simtime; + } looping = false; - dt = target - simtime; } if (predictor and (time1 > 0.0)) { @@ -872,6 +889,7 @@ int SNESSolver::run() { jacobian_pruned = false; // Reset flag. Will be set after pruning. } if (saved_jacobian_lag == 0) { + // This triggers a Jacobian recalculation SNESGetLagJacobian(snes, &saved_jacobian_lag); SNESSetLagJacobian(snes, 1); } @@ -1003,24 +1021,70 @@ int SNESSolver::run() { #endif // PETSC_VERSION_GE(3,20,0) if (looping) { - if (nl_its <= lower_its) { + // Consider changing the timestep. + // Note: The preconditioner depends on the timestep, + // so if it is not recalculated the it will be less + // effective. + if ((nl_its <= lower_its) && (timestep < max_timestep)) { // Increase timestep slightly timestep *= 1.1; if (timestep > max_timestep) { timestep = max_timestep; } + + // Note: Setting the SNESJacobianFn to NULL retains + // previously set evaluation function. + // + // The SNES Jacobian is a combination of the RHS Jacobian + // and a factor involving the timestep. + // Depends on equation_form + // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring); + + // Recompute Jacobian (for now) + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } + } else if (nl_its >= upper_its) { // Reduce timestep slightly timestep *= 0.9; + + // Recompute Jacobian + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } } } } while (looping); + if (!matrix_free) { + ASSERT2(simtime >= target); + ASSERT2(simtime - dt < target); + // Stepped over output timestep => Interpolate + // snes_x is the solution at t = simtime + // x0 is the solution at t = simtime - dt + // Calculate output_x at t = target + VecCopy(snes_x, output_x); + + // Note: If simtime = target then alpha = 0 + // and output_x = snes_x + BoutReal alpha = (simtime - target) / dt; + + // output_x <- alpha * x0 + (1 - alpha) * output_x + VecAXPBY(output_x, alpha, 1. - alpha, x0); + + } else { + // Timestep was adjusted to hit target output time + output_x = snes_x; + } + // Put the result into variables if (scale_vars) { - // scaled_x <- snes_x * var_scaling_factors - int ierr = VecPointwiseMult(scaled_x, snes_x, var_scaling_factors); + // scaled_x <- output_x * var_scaling_factors + int ierr = VecPointwiseMult(scaled_x, output_x, var_scaling_factors); CHKERRQ(ierr); const BoutReal* xdata = nullptr; @@ -1031,15 +1095,15 @@ int SNESSolver::run() { CHKERRQ(ierr); } else { const BoutReal* xdata = nullptr; - int ierr = VecGetArrayRead(snes_x, &xdata); + int ierr = VecGetArrayRead(output_x, &xdata); CHKERRQ(ierr); load_vars(const_cast(xdata)); - ierr = VecRestoreArrayRead(snes_x, &xdata); + ierr = VecRestoreArrayRead(output_x, &xdata); CHKERRQ(ierr); } - run_rhs(simtime); // Run RHS to calculate auxilliary variables + run_rhs(target); // Run RHS to calculate auxilliary variables - if (call_monitors(simtime, s, getNumberOutputSteps()) != 0) { + if (call_monitors(target, s, getNumberOutputSteps()) != 0) { break; // User signalled to quit } } diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 17050ad775..187a64fb40 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -4,7 +4,7 @@ * using PETSc for the SNES interface * ************************************************************************** - * Copyright 2015-2024 BOUT++ contributors + * Copyright 2015-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -116,14 +116,15 @@ private: Vec snes_x; ///< Result of SNES Vec x0; ///< Solution at start of current timestep Vec delta_x; ///< Change in solution + Vec output_x; ///< Solution to output. Used if interpolating. bool predictor; ///< Use linear predictor? Vec x1; ///< Previous solution BoutReal time1{-1.0}; ///< Time of previous solution - SNES snes; ///< SNES context - Mat Jmf; ///< Matrix Free Jacobian - Mat Jfd; ///< Finite Difference Jacobian + SNES snes; ///< SNES context + Mat Jmf; ///< Matrix Free Jacobian + Mat Jfd; ///< Finite Difference Jacobian MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context ///< Jacobian evaluation @@ -135,10 +136,10 @@ private: std::string pc_hypre_type; ///< Hypre preconditioner type std::string line_search_type; ///< Line search type - bool matrix_free; ///< Use matrix free Jacobian - bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? - int lag_jacobian; ///< Re-use Jacobian - bool use_coloring; ///< Use matrix coloring + bool matrix_free; ///< Use matrix free Jacobian + bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? + int lag_jacobian; ///< Re-use Jacobian + bool use_coloring; ///< Use matrix coloring bool jacobian_recalculated; ///< Flag set when Jacobian is recalculated bool prune_jacobian; ///< Remove small elements in the Jacobian? From 8159339f67e2d3f8032ac5ae146e088841879bdf Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 3 May 2025 11:26:30 -0700 Subject: [PATCH 2/2] snes: Make timestep adjustment tunable New tunable inputs for adjustment when timestep fails, and when number of iterations is either too high or too low. --- src/solver/impls/snes/snes.cxx | 33 ++++++++++++++++++++++++--------- src/solver/impls/snes/snes.hxx | 4 ++++ 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 0e514d56cc..780019e5b3 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -123,6 +123,15 @@ SNESSolver::SNESSolver(Options* opts) upper_its((*options)["upper_its"] .doc("Iterations above which the next timestep is reduced") .withDefault(static_cast(maxits * 0.8))), + timestep_factor_on_failure((*options)["timestep_factor_on_failure"] + .doc("Multiply timestep on convergence failure") + .withDefault(0.5)), + timestep_factor_on_upper_its((*options)["timestep_factor_on_upper_its"] + .doc("Multiply timestep if iterations exceed upper_its") + .withDefault(0.9)), + timestep_factor_on_lower_its((*options)["timestep_factor_on_lower_its"] + .doc("Multiply timestep if iterations are below lower_its") + .withDefault(1.4)), diagnose( (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)), diagnose_failures((*options)["diagnose_failures"] @@ -731,6 +740,7 @@ int SNESSolver::run() { bool looping = true; int snes_failures = 0; // Count SNES convergence failures + int steps_since_snes_failure = 0; int saved_jacobian_lag = 0; int loop_count = 0; do { @@ -868,9 +878,10 @@ int SNESSolver::run() { } ++snes_failures; + steps_since_snes_failure = 0; // Try a smaller timestep - timestep /= 2.0; + timestep *= timestep_factor_on_failure; // Restore state VecCopy(x0, snes_x); @@ -961,6 +972,7 @@ int SNESSolver::run() { } simtime += dt; + ++steps_since_snes_failure; if (diagnose) { // Gather and print diagnostic information @@ -970,7 +982,6 @@ int SNESSolver::run() { simtime, timestep, nl_its, lin_its, static_cast(reason)); if (snes_failures > 0) { output.write(", SNES failures: {}", snes_failures); - snes_failures = 0; } output.write("\n"); } @@ -1025,9 +1036,10 @@ int SNESSolver::run() { // Note: The preconditioner depends on the timestep, // so if it is not recalculated the it will be less // effective. - if ((nl_its <= lower_its) && (timestep < max_timestep)) { + if ((nl_its <= lower_its) && (timestep < max_timestep) + && (steps_since_snes_failure > 2)) { // Increase timestep slightly - timestep *= 1.1; + timestep *= timestep_factor_on_lower_its; if (timestep > max_timestep) { timestep = max_timestep; @@ -1041,15 +1053,17 @@ int SNESSolver::run() { // Depends on equation_form // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring); - // Recompute Jacobian (for now) - if (saved_jacobian_lag == 0) { - SNESGetLagJacobian(snes, &saved_jacobian_lag); - SNESSetLagJacobian(snes, 1); + if (static_cast(lin_its) / nl_its > 4) { + // Recompute Jacobian (for now) + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } } } else if (nl_its >= upper_its) { // Reduce timestep slightly - timestep *= 0.9; + timestep *= timestep_factor_on_upper_its; // Recompute Jacobian if (saved_jacobian_lag == 0) { @@ -1058,6 +1072,7 @@ int SNESSolver::run() { } } } + snes_failures = 0; } while (looping); if (!matrix_free) { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 187a64fb40..cba2352582 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -102,6 +102,10 @@ private: int maxits; ///< Maximum nonlinear iterations int lower_its, upper_its; ///< Limits on iterations for timestep adjustment + BoutReal timestep_factor_on_failure; + BoutReal timestep_factor_on_upper_its; + BoutReal timestep_factor_on_lower_its; + bool diagnose; ///< Output additional diagnostics bool diagnose_failures; ///< Print diagnostics on SNES failures