diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 4dce1e7086..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"] @@ -297,13 +306,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 +734,18 @@ 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 steps_since_snes_failure = 0; 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 +823,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)) { @@ -851,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); @@ -872,6 +900,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); } @@ -943,6 +972,7 @@ int SNESSolver::run() { } simtime += dt; + ++steps_since_snes_failure; if (diagnose) { // Gather and print diagnostic information @@ -952,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"); } @@ -1003,24 +1032,74 @@ 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) + && (steps_since_snes_failure > 2)) { // Increase timestep slightly - timestep *= 1.1; + timestep *= timestep_factor_on_lower_its; 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); + + 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) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } } } + snes_failures = 0; } 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 +1110,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..cba2352582 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 * @@ -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 @@ -116,14 +120,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 +140,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?