snes solver: Pseudo Transient Continuation method

[bendudson]

~~Setting pseudo_time = true now enables Pseudo-timestepping, in which each cell has a separate timestep.~~

Setting equation_form = pseudo_transient now enables Pseudo-timestepping, in which each cell has a separate timestep.

The timestep in each cell is set inversely proportional to the residual (RMS time-derivative of all quantities in cell).

Recommend enabling pid_controller that multiplies all timesteps by the same factor, to control the number of nonlinear iterations.

Output time uses the minimum timestep in the domain, so simulations shouldn't need to run for as long as when all cells are evolved at the same speed.

Includes PR #3180

Settings that seem to work well for the Hermes-3 examples/tokamak-1D/1D-threshold case:

[solver]
type = snes
equation_form = pseudo_transient    # Enable pseudo-time stepping

pid_controller = true   # Use PID controller to scale all timesteps based on solver convergence
target_its = 10              # Target number of nonlinear iterations
kP = 0.233
kI = 0.133
kD = 0.0
pid_consider_failures = true  # PID controller becomes cautious after solver failures

snes_type = newtonls
ksp_type = fgmres
pc_type = lu
max_nonlinear_iterations = 20   # Limit significantly above target_its to limit failures
atol = 1e-7
rtol = 1e-6
stol = 1e-12
maxf = 200
maxl = 20

matrix_free_operator = true   # Use the actual nonlinear function for the matrix-vector product

lag_jacobian = 5   # With pid_controller=true the Jacobian will be recalculated every timestep
timestep = 1.0       # Starting internal timestep

[petsc]
snes_ksp_ew = false # Eisenstat-Walker method for setting linear tolerances
ksp_type = fgmres
pc_type = lu
pc_factor_mat_solver_type = strumpack
snes_linesearch_type = basic

This documentation of ADflow contains some useful discussion: https://mdolab-adflow.readthedocs-hosted.com/en/latest/solvers.html The methods used in the ADflow solver seems to be very similar to this: Pseudo-Transient Continuation, matrix-free operator, Inexact Newton-Krylov, Eisenstat-Walker.

When using Eisenstat-Walker the number of nonlinear iterations should be increased: Without it I found target_its=4 was ok, but with EW enabled the solver stalled with small timesteps unless I increased target_its to e.g 10.

Testing Hermes-3 examples/tokamak-1D/1D-threshold on 1 core:

• Previous settings (Hermes-3 master branch) take 1m 18s to run and have not yet reached steady state
• Pseudo-time settings (above) run in 16s and have reached steady state.

For 2D problems try:

[solver]
diagnose = true
type = snes
equation_form = pseudo_transient    # Enable pseudo-time stepping

pid_controller = true   # Use PID controller to scale all timesteps based on solver convergence
target_its = 4         # Target number of nonlinear iterations
kP = 0.233
kI = 0.133
kD = 0.0

snes_type = newtonls
ksp_type = fgmres
pc_type = lu
max_nonlinear_iterations = 20   # Limit significantly above target_its to limit failures
atol = 1e-7
rtol = 1e-6
stol = 1e-12
maxf = 200
maxl = 20

matrix_free_operator = true   # Use the actual nonlinear function for the matrix-vector product

lag_jacobian = 5   # With pid_controller=true the Jacobian will be recalculated every timestep
timestep = 1.0       # Starting internal timestep

[petsc]
snes_ksp_ew = false # Eisenstat-Walker method for setting linear tolerances
ksp_type = fgmres
pc_type = lu
pc_factor_mat_solver_type = strumpack
snes_linesearch_type = basic

[ZedThree]

I've been thinking of expanding the existing solver integrated tests to be able to test different modes of each solver. But because this is really intended for steady-state solutions, we might need a different kind of test -- presumably some function that's constant at large t, and then we only test the final value?

To be clear, the current test does $\int_0^{\pi/2}\sin^2(t) = \pi/4$ (that is, the physics model is $\frac{df}{dt} = \sin^2(t)$ ). The default snes with adaptive = true handles this fine (tolerance of 1e-5). Do we expect this method to be similar?

[bendudson]

Hi @ZedThree ! If there's only one evolving variable then I think this pseudo_transient method should behave like the standard backward Euler method if tolerance is sufficiently small. It's not quite the same though. The real difference comes when there is more than one cell, and each cell starts evolving with its own timestep. I think the better test is something with a known steady state: The time-dependent behavior isn't something we really want to test or constrain, as the sole aim is to get to steady state as quickly as possible.