Using the APT method with auto tuned damping coefficients

This solver is experimental. The DYREL path supports 2D Stokes only. The separate 2D variational-Stokes solver is documented in [2D variational Stokes](variational_stokes.md).

Instead of using the Accelerated Pseudo-Transient method, where the damping coefficients are constant throughout the PT iterations (as in Räss et al., 2022), we can use a self-tuning version of the APT method based on the approach described in Duretz et al., 2025.

Usage

To use this solver, only two changes are needed with respect to the scripts using the APT solver described in previous examples:

1. The PTStokesCoeffs object containing the arrays needed for the standard APT solver is not needed anymore, and needs to be replaced by the DYREL object that contains all the new arrays that are needed for the self-tuned APT method. This means we need to change this line

pt_stokes = PTStokesCoeffs(li, di; ϵ_abs = 1.0e-6, ϵ_rel = 1.0e-6, CFL = 0.95 / √2)

by this one

dyrel = DYREL(backend, stokes, rheology, phase_ratios, di, dt; ϵ=1e-6)

Note that the `DYREL` arrays need the effective viscosity of the model, so `DYREL` must be instantiated _after_ an effective-viscosity guess is available.

1. The last change requires changing the solver function call to the following:

solve_DYREL!(
    stokes,
    ρg,
    dyrel,
    flow_bcs,
    phase_ratios,
    rheology,
    args,
    di,
    dt,
    igg;
    kwargs = (;
        iterMax              = 50.0e3,
        nout                 = 10,
        rel_drop             = 0.1,
        λ_relaxation_PH      = 1,
        λ_relaxation_DR      = 1,
        verbose_PH           = false,
        verbose_DR           = false,
        viscosity_relaxation = 1,
        linear_viscosity     = true,
        viscosity_cutoff     = (-Inf, Inf),
        free_surface         = false,
    )
);

where the solver keyword arguments are:

• iterMax $\to$ maximum number of total iterations.

• nout $\to$ damping coefficients are re-computed every nout iterations.

• rel_drop $\to$ the tolerance for the inner dynamic relaxation loop is $\text{error(P^n) text{rel\_drop}}$ where $n$ is the inner Powell-Hesteness iteration counter.

• λ_relaxation_PH $\to$ relaxation coefficient for the plastic multiplier ($\cdot \lambda$) during the inner Powell-Hesteness loop. λ_relaxation_PH=1 means no relaxation.

• λ_relaxation_DR $\to$ relaxation coefficient for the plastic multiplier ($\cdot \lambda$) during the inner Dynamic Relaxation loop. λ_relaxation_DR=1 means no relaxation.

• verbose_PH $\to$ print solver metrics during the inner Powell–Hestenes loop.

• verbose_DR $\to$ print solver metrics during the inner Dynamic Relaxation loop.

• viscosity_relaxation $\to$ relaxation coefficient for the viscosity. viscosity_relaxation=1 means no relaxation.

• linear_viscosity $\to$ if the rheology is linear (viscosity will not be updated during the solver iterations).

• viscosity_cutoff $\to$ viscosity is clamped so that $\text{text{viscosity\_cutoff}\_1 leq eta leq text{viscosity\_cutoff}\_2}$.

• free_surface $\to$ include the 2D density-gradient free-surface stabilization in both the momentum residual and the self-tuned vertical pseudo-transient coefficients. The default is false.

When free_surface=true, DYREL adds the local diagonal $-\mathrm{\Delta }t{\partial }_y(\rho g_y)$ to Dy and to the corresponding Gershgorin row bound whenever the pseudo-transient coefficients are refreshed. The same term is used by the Powell–Hestenes and dynamic-relaxation residual kernels.

2D variational DYREL

Free surfaces represented by a RockRatio use the dedicated 2D entry point solve_VariationalDYREL!. It combines DYREL's Powell–Hestenes and dynamic relaxation iteration with the volume-weighted operator of Larionov, Batty, and Bridson (2017): center fractions weight pressure and normal stress, vertex fractions weight shear stress, and face fractions weight momentum rows.

ϕ = RockRatio(backend, ni)
update_rock_ratio!(ϕ, phase_ratios, air_phase)
dyrel = DYREL(backend, stokes, rheology, phase_ratios, ϕ, grid.di, dt; ϵ = 1.0e-6)

solve_VariationalDYREL!(
    stokes, ρg, dyrel, flow_bcs, phase_ratios, ϕ,
    rheology, args, grid, dt, igg;
    air_phase = air_phase,
    linear_viscosity = true,
    free_surface = true,
    pressure_relaxation = 0.5,
    verbose_PH = false,
    verbose_DR = false,
)

The constructor and solver must receive the same RockRatio. Zero-volume pressure and velocity rows are eliminated, positive sliver faces use a bounded face mass, and changing the mask between calls resets the dynamic-relaxation history. air_phase excludes the air phase from material averages, while pressure_relaxation damps the Powell–Hestenes pressure update; its default is 1. The solver accepts either a Geometry{2} or a legacy uniform-spacing tuple, and takes its options either as the plain keywords above or bundled as a single kwargs = (; ...) NamedTuple. Standard solve_DYREL! remains unchanged for unweighted 2D problems.

A velocity row whose preconditioner diagonal comes out zero is decoupled from its neighbours and is preconditioned with the identity. A diagonal that comes out negative or non-finite is not: it means the free-surface term is large enough to invert the row, or that a degenerate phase sample produced a non-finite viscosity, and it is propagated as NaN so the solve reports a failure rather than continuing with an arbitrary substitute. Reaching it with free_surface = true indicates the timestep is too large for the density contrast at the surface.

See 2D variational Stokes for marker-chain ordering, mask construction, and the mathematical reference.

Examples

Examples of a set of miniapps using this solver can be found in this folder.