Skip to content

Use SNESSetLagJacobianPersists to force recalculating Jacobian #3404

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cmacmackin wants to merge 1 commit into
base: next
Choose a base branch
from
Open

Use SNESSetLagJacobianPersists to force recalculating Jacobian #3404

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 21 additions & 31 deletions src/solver/impls/snes/snes.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -881,7 +881,7 @@ int SNESSolver::init() {
return 0;
}

PetscErrorCode SNESSolver::rescale(int& saved_jacobian_lag) {
PetscErrorCode SNESSolver::rescale() {
// Individual variable scaling
// Note: If variables are rescaled then the Jacobian columns
// need to be scaled or recalculated
Expand Down Expand Up @@ -929,12 +929,8 @@ 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);
PetscCall(SNESSetLagJacobianPersists(snes, PETSC_FALSE));

return PETSC_SUCCESS;
}

Expand All @@ -957,9 +953,8 @@ int SNESSolver::run() {
PetscCall(VecRestoreArray(snes_x, &xdata));
}

int saved_jacobian_lag = 0;
if (scale_vars) {
PetscCall(rescale(saved_jacobian_lag));
PetscCall(rescale());
}

// Initialise residuals
Expand All @@ -980,6 +975,8 @@ int SNESSolver::run() {
recent_failure_rate = 0.0;
int steps_since_rescale = 0;
BoutReal change_since_rescale = 0.;
// Tracks whether we force persisting of lagged Jacobian to be turned off for the next time-step
bool persist_jac_lag_off = false;
for (int s = 1; s <= getNumberOutputSteps(); s++) {
target += getOutputTimestep();

Expand All @@ -994,10 +991,17 @@ int SNESSolver::run() {
break; // Could happen if step over multiple outputs
}

if (persist_jac_lag_off) {
persist_jac_lag_off = false;
} else {
PetscCall(
SNESSetLagJacobianPersists(snes, static_cast<PetscBool>(jacobian_persists)));
}

if (scale_vars
and (change_since_rescale > rescale_threshold
or steps_since_rescale == rescale_period)) {
PetscCall(rescale(saved_jacobian_lag));
PetscCall(rescale());
change_since_rescale = 0.;
steps_since_rescale = 0;
}
Expand Down Expand Up @@ -1026,14 +1030,14 @@ int SNESSolver::run() {
if (timestep < dt_min_reset) {
// Hit the minimum timestep, probably through repeated failures

if (saved_jacobian_lag != 0) {
if (persist_jac_lag_off) {
// Already tried this and it didn't work
throw BoutException("Solver failed after many attempts");
}

// Try resetting the preconditioner, turn off predictor, and use a large timestep
SNESGetLagJacobian(snes, &saved_jacobian_lag);
SNESSetLagJacobian(snes, 1);
PetscCall(SNESSetLagJacobianPersists(snes, PETSC_FALSE));
persist_jac_lag_off = true;
timestep = getOutputTimestep();
predictor = false; // Predictor can cause problems in near steady-state.
}
Expand All @@ -1060,8 +1064,6 @@ int SNESSolver::run() {
const BoutReal beta = dt / (simtime - time1);
VecAXPBY(snes_x, -beta, (1. + beta), x1);
}

SNESSetLagJacobian(snes, lag_jacobian);
}

// Run the solver
Expand Down Expand Up @@ -1142,11 +1144,8 @@ int SNESSolver::run() {
PetscCall(FDJrestoreFromPruning());
}

if (saved_jacobian_lag == 0) {
// This triggers a Jacobian recalculation
SNESGetLagJacobian(snes, &saved_jacobian_lag);
SNESSetLagJacobian(snes, 1);
}
PetscCall(SNESSetLagJacobianPersists(snes, PETSC_FALSE));
persist_jac_lag_off = true;

// Check lock state
PetscInt lock_state;
Expand All @@ -1162,13 +1161,6 @@ int SNESSolver::run() {
continue; // Try again
}

if (saved_jacobian_lag != 0) {

@bendudson bendudson Jun 23, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should restore the previous Jacobian lag on successful SNES convergence. I think this should restore the JacobianPersists setting to its input option. It looks like rescale() will set JacobianPersists to False, but where would it be restored to True (if that's what was in the settings)?

@cmacmackin cmacmackin Jun 24, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If rescale determines that JacobianPersists should be false then that applies for the current time-step and it will be reset on the next one.

@bendudson bendudson Jun 24, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where does the reset happen? This point (around line 1165) is where the Jacobian lag was being reset: Execution reaches this point only after a successful SNES convergence. If this is removed then I don't see where JacobianPersists is reset.

// Following successful step, reset Jacobian lag
// to previous value
SNESSetLagJacobian(snes, saved_jacobian_lag);
saved_jacobian_lag = 0;
}

if (predictor) {
// Save previous values: x1 <- x0
VecCopy(x0, x1);
Expand Down Expand Up @@ -1275,10 +1267,8 @@ int SNESSolver::run() {

if (static_cast<BoutReal>(lin_its) / nl_its > 0.5 * maxl) {
// Recompute Jacobian if number of linear iterations is too high
if (saved_jacobian_lag == 0) {
SNESGetLagJacobian(snes, &saved_jacobian_lag);
SNESSetLagJacobian(snes, 1);
}
PetscCall(SNESSetLagJacobianPersists(snes, PETSC_FALSE));
persist_jac_lag_off = true;
}

snes_failures = 0;
Expand Down
2 changes: 1 addition & 1 deletion src/solver/impls/snes/snes.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ private:

/// 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);
PetscErrorCode rescale();

/// Call the physics model RHS function
///
Expand Down
Loading