Skip to content

snes: pseudo_alpha_minimum setting #3377

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
bendudson wants to merge 1 commit into
base: next
Choose a base branch
from
Open

snes: pseudo_alpha_minimum setting #3377

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
42 changes: 27 additions & 15 deletions src/solver/impls/snes/snes.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,18 +4,25 @@

#include "snes.hxx"

#include <bout/assert.hxx>
#include <bout/bout_types.hxx>
#include <bout/boutcomm.hxx>
#include <bout/boutexception.hxx>
#include <bout/field.hxx>
#include <bout/field3d.hxx>
#include <bout/globals.hxx>
#include <bout/output.hxx>
#include <bout/output_bout_types.hxx>
#include <bout/petsc_interface.hxx>
#include <bout/solver.hxx>
#include <bout/utils.hxx>
#include <bout/unused.hxx>

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <set>
#include <utility>
#include <vector>

#include "petscerror.h"
Expand Down Expand Up @@ -306,7 +313,7 @@ PetscErrorCode SNESSolver::FDJinitialise() {
}
// 3D fields
for (int z = mesh->zstart; z <= mesh->zend; z++) {
int ind = ROUND(index(x, y, z));
const int ind = ROUND(index(x, y, z));

for (int i = 0; i < n3d; i++) {
PetscInt row = ind + i;
Expand Down Expand Up @@ -342,7 +349,7 @@ PetscErrorCode SNESSolver::FDJinitialise() {
// 3D fields on this cell
for (int j = 0; j < n3d; j++) {
const PetscInt col = ind2 + j;
PetscErrorCode ierr =
const PetscErrorCode ierr =
MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES);

if (ierr != PETSC_SUCCESS) {
Expand Down Expand Up @@ -536,6 +543,9 @@ SNESSolver::SNESSolver(Options* opts)
pseudo_alpha((*options)["pseudo_alpha"]
.doc("Sets timestep using dt = alpha / residual")
.withDefault(100. * atol * timestep)),
pseudo_alpha_minimum(
(*options)["pseudo_alpha_minimum"].doc("Minimum value of pseudo_alpha")
.withDefault(0.1 * pseudo_alpha)),
pseudo_growth_factor((*options)["pseudo_growth_factor"]
.doc("PTC growth factor on success")
.withDefault(1.1)),
Expand Down Expand Up @@ -884,7 +894,7 @@ int SNESSolver::run() {
int saved_jacobian_lag = 0;
int loop_count = 0;

BoutReal start_global_residual = global_residual;
const BoutReal start_global_residual = global_residual;
do {
if ((output_trigger == BoutSnesOutput::fixed_time_interval && (simtime >= target))
|| (output_trigger == BoutSnesOutput::residual_ratio
Expand Down Expand Up @@ -1181,8 +1191,10 @@ int SNESSolver::run() {

if (equation_form == BoutSnesEquationForm::pseudo_transient) {
// Adjust pseudo_alpha to globally scale timesteps
pseudo_alpha = updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate,
max_timestep * atol * 100);
pseudo_alpha = std::max({
updateGlobalTimestep(pseudo_alpha, nl_its, recent_failure_rate,
max_timestep * atol * 100),
pseudo_alpha_minimum});

// Adjust local timesteps
PetscCall(updatePseudoTimestepping());
Expand Down Expand Up @@ -1216,7 +1228,7 @@ int SNESSolver::run() {

// Note: If simtime = target then alpha = 0
// and output_x = snes_x
BoutReal alpha = (simtime - target) / dt;
const BoutReal alpha = (simtime - target) / dt;

// output_x <- alpha * x0 + (1 - alpha) * output_x
VecAXPBY(output_x, alpha, 1. - alpha, x0);
Expand Down Expand Up @@ -1437,7 +1449,7 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() {
}
if (count > 0) {
// Adjust timestep for these quantities
BoutReal new_timestep = updatePseudoTimestep(
const BoutReal new_timestep = updatePseudoTimestep(
dt_data[idx], local_residual_2d_prev[i2d], local_residual_2d[i2d]);
for (int i = 0; i != count; ++i) {
dt_data[idx++] = new_timestep;
Expand Down Expand Up @@ -1529,9 +1541,10 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping() {
/// rapid changes in timestep.
BoutReal SNESSolver::updatePseudoTimestep_inverse_residual(BoutReal previous_timestep,
BoutReal current_residual) {
return std::min(
{std::max({pseudo_alpha / current_residual, previous_timestep / 1.5, dt_min_reset}),
1.5 * previous_timestep, max_timestep});
return std::max({
std::min({std::max({pseudo_alpha / current_residual, previous_timestep / 1.5}),
1.5 * previous_timestep, max_timestep}),
dt_min_reset});
}

// Strategy based on history of residuals
Expand All @@ -1552,7 +1565,8 @@ BoutReal SNESSolver::updatePseudoTimestep_history_based(BoutReal previous_timest
if (reduction_ratio < 0.8) {
return std::min(0.5 * (pseudo_growth_factor + 1.) * previous_timestep,
max_timestep);
} else if (reduction_ratio > 1.2) {
}
if (reduction_ratio > 1.2) {
return std::max(0.5 * (pseudo_reduction_factor + 1) * previous_timestep,
dt_min_reset);
}
Expand Down Expand Up @@ -1780,10 +1794,8 @@ PetscErrorCode SNESSolver::scaleJacobian(Mat Jac_new) {
// Calculate a norm of this row of the Jacobian
PetscScalar norm = 0.0;
for (int col = 0; col < ncols; col++) {
PetscScalar absval = std::abs(vals[col]);
if (absval > norm) {
norm = absval;
}
const PetscScalar absval = std::abs(vals[col]);
norm = std::max({norm, absval});
// Can we identify small elements and remove them?
// so we don't need to calculate them next time
}
Expand Down
1 change: 1 addition & 0 deletions src/solver/impls/snes/snes.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,6 +153,7 @@ private:
// These are used if equation_form = pseudo_transient
BoutPTCStrategy pseudo_strategy; ///< Strategy to use when setting timesteps
BoutReal pseudo_alpha; ///< dt = alpha / residual
BoutReal pseudo_alpha_minimum; ///< Minimum value of alpha
BoutReal pseudo_growth_factor; ///< Timestep increase 1.1 - 1.2
BoutReal pseudo_reduction_factor; ///< Timestep decrease 0.5
BoutReal pseudo_max_ratio; ///< Maximum timestep ratio between neighboring cells
Expand Down
Loading