11.7. The SUNNonlinSol_Newton implementation

This section describes the SUNNonlinSol implementation of Newton’s method. To access the SUNNonlinSol_Newton module, include the header file sunnonlinsol/sunnonlinsol_newton.h. We note that the SUNNonlinSol_Newton module is accessible from SUNDIALS integrators without separately linking to the libsundials_sunnonlinsolnewton module library.

11.7.1. SUNNonlinSol_Newton description

To find the solution to

(11.9)\[F(y) = 0\]

given an initial guess \(y^{(0)}\), Newton’s method computes a series of approximate solutions

\[y^{(m+1)} = y^{(m)} + \delta^{(m+1)}\]

where \(m\) is the Newton iteration index, and the Newton update \(\delta^{(m+1)}\) is the solution of the linear system

(11.10)\[A(y^{(m)}) \delta^{(m+1)} = -F(y^{(m)}) \, ,\]

in which \(A\) is the Jacobian matrix

(11.11)\[A \equiv \partial F / \partial y \, .\]

Depending on the linear solver used, the SUNNonlinSol_Newton module will employ either a Modified Newton method or an Inexact Newton method [24, 29, 46, 48, 95]. When used with a direct linear solver, the Jacobian matrix \(A\) is held constant during the Newton iteration, resulting in a Modified Newton method. With a matrix-free iterative linear solver, the iteration is an Inexact Newton method.

In both cases, calls to the integrator-supplied SUNNonlinSolLSetupFn function are made infrequently to amortize the increased cost of matrix operations (updating \(A\) and its factorization within direct linear solvers, or updating the preconditioner within iterative linear solvers). Specifically, SUNNonlinSol_Newton will call the SUNNonlinSolLSetupFn function in two instances:

when requested by the integrator (the input callLSetSetup is SUNTRUE) before attempting the Newton iteration, or
when reattempting the nonlinear solve after a recoverable failure occurs in the Newton iteration with stale Jacobian information (jcur is SUNFALSE). In this case, SUNNonlinSol_Newton will set jbad to SUNTRUE before calling the SUNNonlinSolLSetupFn() function.

Whether the Jacobian matrix \(A\) is fully or partially updated depends on logic unique to each integrator-supplied SUNNonlinSolLSetupFn routine. We refer to the discussion of nonlinear solver strategies provided in the package-specific Mathematics section of the documentation for details.

The default maximum number of iterations and the stopping criteria for the Newton iteration are supplied by the SUNDIALS integrator when SUNNonlinSol_Newton is attached to it. Both the maximum number of iterations and the convergence test function may be modified by the user by calling the SUNNonlinSolSetMaxIters() and/or SUNNonlinSolSetConvTestFn() functions after attaching the SUNNonlinSol_Newton object to the integrator.

11.7.2. SUNNonlinSol_Newton functions

The SUNNonlinSol_Newton module provides the following constructor for creating the SUNNonlinearSolver object.

SUNNonlinearSolver SUNNonlinSol_Newton(N_Vector y, SUNContext sunctx)

This creates a SUNNonlinearSolver object for use with SUNDIALS integrators to solve nonlinear systems of the form \(F(y) = 0\) using Newton’s method.

Arguments:

y – a template for cloning vectors needed within the solver.
sunctx – the SUNContext object (see §1.3)

Return value:

A SUNNonlinSol object if the constructor exits successfully, otherwise it will be NULL.

The SUNNonlinSol_Newton module implements all of the functions defined in §11.1.1–§11.1.3 except for SUNNonlinSolSetup(). The SUNNonlinSol_Newton functions have the same names as those defined by the generic SUNNonlinSol API with _Newton appended to the function name. Unless using the SUNNonlinSol_Newton module as a standalone nonlinear solver the generic functions defined in §11.1.1–§11.1.3 should be called in favor of the SUNNonlinSol_Newton-specific implementations.

The SUNNonlinSol_Newton module also defines the following user-callable functions.

SUNErrCode SUNNonlinSolGetSysFn_Newton(SUNNonlinearSolver NLS, SUNNonlinSolSysFn *SysFn)

This returns the residual function that defines the nonlinear system.

Arguments:

NLS – a SUNNonlinSol object.
SysFn – the function defining the nonlinear system.

Return value:

A SUNErrCode

Notes:

This function is intended for users that wish to evaluate the nonlinear residual in a custom convergence test function for the SUNNonlinSol_Newton module. We note that SUNNonlinSol_Newton will not leverage the results from any user calls to SysFn.

SUNErrCode SUNNonlinSolSetComputeStiffnessRatio_Newton(SUNNonlinearSolver NLS, sunbooleantype onoff)

This enables or disables the additional residual norm evaluation used to compute the stiffness ratio stiffr from [113].

Arguments:

NLS – a SUNNonlinSol object.
onoff – SUNTRUE to compute stiffr after each continued Newton iteration, or SUNFALSE to disable that extra work.

Return value:

A SUNErrCode

Notes:

By default the stiffness ratio computation is disabled. This is automatically enabled when the Newton solver is used from within the SUNNonlinSol_Auto module (see §11.9).

Added in version 7.8.0.

SUNErrCode SUNNonlinSolGetStiffnessRatio_Newton(SUNNonlinearSolver NLS, sunrealtype *stiffr)

This returns the most recently computed stiffness ratio.

Arguments:

NLS – a SUNNonlinSol object.
stiffr – the stiffness metric \(\|F(y^m)\| / \|\delta_m\|\).

Return value:

A SUNErrCode

Added in version 7.8.0.

11.7.3. SUNNonlinSol_Newton content

The content field of the SUNNonlinSol_Newton module is the following structure.

struct _SUNNonlinearSolverContent_Newton {

  SUNNonlinSolSysFn      Sys;
  SUNNonlinSolLSetupFn   LSetup;
  SUNNonlinSolLSolveFn   LSolve;
  SUNNonlinSolConvTestFn CTest;
  SUNNonlinSolNormFn     norm_fn;
  void                  *norm_fn_data;

  N_Vector       delta;
  sunbooleantype jcur;
  int            curiter;
  int            maxiters;
  long int       niters;
  long int       nconvfails;
  sunbooleantype compute_stiffr;
  sunrealtype    stiffr;
  sunrealtype    delnrm;
  void*          ctest_data;
};

These entries of the content field contain the following information:

Sys – the function for evaluating the nonlinear system,
LSetup – the package-supplied function for setting up the linear solver,
LSolve – the package-supplied function for performing a linear solve,
CTest – the function for checking convergence of the Newton iteration,
norm_fn – an optional callback for computing the update norm or residual norm,
norm_fn_data – the data pointer passed to norm_fn,
delta – the Newton iteration update vector,
jcur – the Jacobian status (SUNTRUE = current, SUNFALSE = stale),
curiter – the current number of iterations in the solve attempt,
maxiters – the maximum number of Newton iterations allowed in a solve,
niters – the total number of nonlinear iterations across all solves,
nconvfails – the total number of nonlinear convergence failures across all solves,
compute_stiffr – a flag indicating whether to compute the stiffness ratio after continued iterations,
stiffr – the most recently computed stiffness ratio,
delnrm – the WRMS norm of the most recent Newton update,
ctest_data – the data pointer passed to the convergence test function,

11.8. The SUNNonlinSol_FixedPoint implementation

This section describes the SUNNonlinSol implementation of a fixed point (functional) iteration with optional Anderson acceleration. To access the SUNNonlinSol_FixedPoint module, include the header file sunnonlinsol/sunnonlinsol_fixedpoint.h. We note that the SUNNonlinSol_FixedPoint module is accessible from SUNDIALS integrators without separately linking to the libsundials_sunnonlinsolfixedpoint module library.

11.8.1. SUNNonlinSol_FixedPoint description

To find the solution to

(11.12)\[G(y) = y \,\]

given an initial guess \(y^{(0)}\), the fixed point iteration computes a series of approximate solutions

(11.13)\[y^{(n+1)} = G(y^{(n)})\]

where \(n\) is the iteration index. The convergence of this iteration may be accelerated using Anderson’s method [11, 58, 108, 166]. With Anderson acceleration using subspace size \(m\), the series of approximate solutions can be formulated as the linear combination

(11.14)\[y^{(n+1)} = \beta \sum_{i=0}^{m_n} \alpha_i^{(n)} G(y^{(n-m_n+i)}) + (1 - \beta) \sum_{i=0}^{m_n} \alpha_i^{(n)} y_{n-m_n+i}\]

where \(m_n = \min{\{m,n\}}\) and the factors

\[\alpha^{(n)} =(\alpha_0^{(n)}, \ldots, \alpha_{m_n}^{(n)})\]

solve the minimization problem \(\min\limits_\alpha \| F_n \alpha^T \|_2\) under the constraint that \(\sum\limits_{i=0}^{m_n} \alpha_i = 1\) where

\[F_{n} = (f_{n-m_n}, \ldots, f_{n})\]

with \(f_i = G(y^{(i)}) - y^{(i)}\). Due to this constraint, in the limit of \(m=0\) the accelerated fixed point iteration formula (11.14) simplifies to the standard fixed point iteration (11.13).

Following the recommendations made in [166], the SUNNonlinSol_FixedPoint implementation computes the series of approximate solutions as

(11.15)\[y^{(n+1)} = G(y^{(n)})-\sum_{i=0}^{m_n-1} \gamma_i^{(n)} \Delta g_{n-m_n+i} - (1 - \beta) (f(y^{(n)}) - \sum_{i=0}^{m_n-1} \gamma_i^{(n)} \Delta f_{n-m_n+i})\]

with \(\Delta g_i = G(y^{(i+1)}) - G(y^{(i)})\) and where the factors

\[\gamma^{(n)} =(\gamma_0^{(n)}, \ldots, \gamma_{m_n-1}^{(n)})\]

solve the unconstrained minimization problem \(\min\limits_\gamma \| f_n - \Delta F_n \gamma^T \|_2\) where

\[\Delta F_{n} = (\Delta f_{n-m_n}, \ldots, \Delta f_{n-1}),\]

with \(\Delta f_i = f_{i+1} - f_i\). The least-squares problem is solved by applying a QR factorization to \(\Delta F_n = Q_n R_n\) and solving \(R_n \gamma = Q_n^T f_n\).

The acceleration subspace size \(m\) is required when constructing the SUNNonlinSol_FixedPoint object. The default maximum number of iterations and the stopping criteria for the fixed point iteration are supplied by the SUNDIALS integrator when SUNNonlinSol_FixedPoint is attached to it. Both the maximum number of iterations and the convergence test function may be modified by the user by calling SUNNonlinSolSetMaxIters() and SUNNonlinSolSetConvTestFn() after attaching the SUNNonlinSol_FixedPoint object to the integrator.

11.8.2. SUNNonlinSol_FixedPoint functions

The SUNNonlinSol_FixedPoint module provides the following constructor for creating the SUNNonlinearSolver object.

SUNNonlinearSolver SUNNonlinSol_FixedPoint(N_Vector y, int m, SUNContext sunctx)

This creates a SUNNonlinearSolver object for use with SUNDIALS integrators to solve nonlinear systems of the form \(G(y) = y\).

Arguments:

y – a template for cloning vectors needed within the solver.
m – the number of acceleration vectors to use.
sunctx – the SUNContext object (see §1.3)

Return value:

A SUNNonlinSol object if the constructor exits successfully, otherwise it will be NULL.

Since the accelerated fixed point iteration (11.13) does not require the setup or solution of any linear systems, the SUNNonlinSol_FixedPoint module implements all of the functions defined in §11.1.1–§11.1.3 except for the SUNNonlinSolSetup(), SUNNonlinSolSetLSetupFn(), and SUNNonlinSolSetLSolveFn() functions, that are set to NULL. The SUNNonlinSol_FixedPoint functions have the same names as those defined by the generic SUNNonlinSol API with _FixedPoint appended to the function name. Unless using the SUNNonlinSol_FixedPoint module as a standalone nonlinear solver the generic functions defined in §11.1.1–§11.1.3 should be called in favor of the SUNNonlinSol_FixedPoint-specific implementations.

The SUNNonlinSol_FixedPoint module also defines the following user-callable functions.

SUNErrCode SUNNonlinSolGetSysFn_FixedPoint(SUNNonlinearSolver NLS, SUNNonlinSolSysFn *SysFn)

This returns the fixed-point function that defines the nonlinear system.

Arguments:

NLS – a SUNNonlinSol object.
SysFn – the function defining the nonlinear system.

Return value:

A SUNErrCode

Notes:

This function is intended for users that wish to evaluate the fixed-point function in a custom convergence test function for the SUNNonlinSol_FixedPoint module. We note that SUNNonlinSol_FixedPoint will not leverage the results from any user calls to SysFn.

SUNErrCode SUNNonlinSolSetDamping_FixedPoint(SUNNonlinearSolver NLS, sunrealtype beta)

This sets the damping parameter \(\beta\) to use with Anderson acceleration. By default damping is disabled i.e., \(\beta = 1.0\).

Arguments:

NLS – a SUNNonlinSol object.
beta – the damping parameter \(0 < \beta \leq 1\).

Return value:

A SUNErrCode

Notes:

A beta value should satisfy \(0 < \beta < 1\) if damping is to be used. A value of one or more will disable damping.

If supported by the SUNNonlinearSolver implementation, this routine will be called by SUNNonlinSolSetOptions() when using the key NLSid.damping.

11.8.3. SUNNonlinSol_FixedPoint content

The content field of the SUNNonlinSol_FixedPoint module is the following structure.

struct _SUNNonlinearSolverContent_FixedPoint {

  SUNNonlinSolSysFn      Sys;
  SUNNonlinSolConvTestFn CTest;
  SUNNonlinSolNormFn     norm_fn;
  void                  *norm_fn_data;

  int            m;
  int           *imap;
  sunrealtype   *R;
  sunbooleantype damping;
  sunrealtype    beta;
  sunrealtype    *gamma;
  sunrealtype    *cvals;
  sunrealtype    delnrm;
  sunrealtype    crate;
  sunrealtype    crate_const;
  N_Vector      *df;
  N_Vector      *dg;
  N_Vector      *q;
  N_Vector      *Xvecs;
  N_Vector       yprev;
  N_Vector       gy;
  N_Vector       fold;
  N_Vector       gold;
  N_Vector       delta;
  int            curiter;
  int            maxiters;
  long int       niters;
  long int       nconvfails;
  void          *ctest_data;
};

The following entries of the content field are always allocated:

Sys – function for evaluating the nonlinear system,
CTest – function for checking convergence of the fixed point iteration,
norm_fn – optional callback for computing the update norm or residual norm,
norm_fn_data – user data passed to norm_fn,
yprev – N_Vector used to store previous fixed-point iterate,
gy – N_Vector used to store \(G(y)\) in fixed-point algorithm,
delta – N_Vector used to store difference between successive fixed-point iterates,
delnrm – WRMS norm of the most recent fixed-point update,
crate – estimated nonlinear convergence rate,
crate_const – convergence rate constant used in the crate estimate,
curiter – the current number of iterations in the solve attempt,
maxiters – the maximum number of fixed-point iterations allowed in a solve,
niters – the total number of nonlinear iterations across all solves,
nconvfails – the total number of nonlinear convergence failures across all solves,
ctest_data – the data pointer passed to the convergence test function,
m – number of acceleration vectors,

If Anderson acceleration is requested (i.e., \(m>0\) in the call to SUNNonlinSol_FixedPoint()), then the following items are also allocated within the content field:

imap – index array used in acceleration algorithm (length m),
damping – a flag indicating if damping is enabled,
beta – the damping parameter,
R – small matrix used in acceleration algorithm (length m*m),
gamma – small vector used in acceleration algorithm (length m),
cvals – small vector used in acceleration algorithm (length m+1),
df – array of vectors used in acceleration algorithm (length m),
dg – array of vectors used in acceleration algorithm (length m),
q – array of vectors used in acceleration algorithm (length m),
Xvecs – vector pointer array used in acceleration algorithm (length m+1),
fold – vector used in acceleration algorithm, and
gold – vector used in acceleration algorithm.

11.9. The SUNNonlinSol_Auto implementation

Added in version 7.8.0.

This section describes the SUNNonlinSol implementation that can automatically switch between §11.8 and §11.7 during a solve. The switching algorithm is based on [113].

To access the SUNNonlinSol_Auto module, include the header file sunnonlinsol/sunnonlinsol_auto.h. The library to link to is libsundials_sunnonlinsolauto.lib where .lib is typically .so for shared libraries and .a for static libraries.

11.9.1. SUNNonlinSol_Auto description

SUNNonlinSol_Auto is a hybrid nonlinear solver that delegates each nonlinear solve to an underlying fixed-point or Newton solver and may request a switch to the other algorithm based on a runtime switching criterion.

Switching decisions are checked at every nonlinear iteration by intercepting the convergence test function provided by the integrator (see SUNNonlinSolSetConvTestFn()). When SUNNonlinSol_Auto decides to switch, the active sub-solver returns SUN_NLS_SWITCH from the convergence test. SUNNonlinSol_Auto handles this internally by retrying the current nonlinear solve with the other sub-solver, so the integrator does not reduce the step size solely because the nonlinear algorithm changed.

A full mathematical description of the switching criterion and algorithm can be found in [113]. In short, switching from fixed-point to Newton occurs when the fixed-point convergence-rate estimate supplied through SUNNonlinSolSetGetConvRateFn()

\[R \leftarrow \max\{0.3R, \|\delta_m\| / \|\delta_{m-1}\|\},\]

indicates slow convergence or divergence, i.e., when \(R \ge \alpha\), for a specified number of nonlinear solves. The implementation default is \(\alpha = 0.8\) with fp_to_newt_delay = 1, consistent with the Norsett switching criterion. Switching from Newton to fixed-point occurs when the stiffness indicator

\[\text{stiffr} \leftarrow \|F(y^m)\| / \|\delta_m\|,\]

satisfies \(\text{stiffr} < \beta\) for a specified number of nonlinear solves. The implementation default is \(\beta = 2.0\) with newt_to_fp_delay = 10 solves before switching from Newton to fixed-point is allowed, matching the recommendation in [113].

11.9.2. SUNNonlinSol_Auto functions

The SUNNonlinSol_Auto module provides the following constructor for creating the SUNNonlinearSolver object.

SUNNonlinearSolver SUNNonlinSol_Auto(N_Vector y, int m, SUNNonlinSolAutoType initial_solver_type, SUNContext sunctx)

This creates a SUNNonlinearSolver object for use with SUNDIALS integrators.

Parameters:

y – a template for cloning vectors needed within the solver.
m – the number of acceleration vectors to use with the underlying fixed-point solver (passed to SUNNonlinSol_FixedPoint()).
initial_solver_type – the initial solver type.
sunctx – the SUNContext object (see §1.3)

Returns:

a pointer to a SUNNonlinearSolver object if the constructor exits successfully, otherwise it will be NULL.

enum SUNNonlinSolAutoType

An identifier indicating which underlying solver is used initially.

enumerator SUNNONLINSOL_AUTO_FIXEDPOINT: Use the fixed-point solver.

enumerator SUNNONLINSOL_AUTO_NEWTON: Use the Newton solver.

The SUNNonlinSol_Auto module follows the generic nonlinear-solver interface defined in §11.1.1–§11.1.3. Users should normally call the generic SUNNonlinSol API. When solver-specific access is needed, the auto module provides _Auto helper routines for switching parameters, accumulated statistics, the currently active solver type, and access to the underlying fixed-point and Newton solvers. The generic SUNNonlinSolSetMaxIters() routine forwards the same nonlinear-iteration limit to both wrapped sub-solvers. Likewise, the generic SUNNonlinSolGetNumConvFails() routine returns the sum of the wrapped solvers’ convergence failures from the most recent Auto solve. The generic SUNNonlinSolSetOptions() routine only applies the Auto-specific switching_parameters key and the generic max_iters key directly to the Auto solver; the max_iters key therefore updates both wrapped sub-solvers. To configure other fixed-point- or Newton-specific options separately, retrieve the wrapped solver with the getter functions below and call the appropriate module-specific setter on that object. Likewise, calling the generic SUNNonlinSolSetNormFn(), SUNNonlinSolSetGetUpdateNormFn(), or SUNNonlinSolSetGetConvRateFn() on the Auto solver forwards that callback to the wrapped sub-solvers so they use the same integrator-provided convergence information.

The SUNNonlinSol_Auto module also defines the following function for controlling the switching behavior.

SUNErrCode SUNNonlinSolSetSwitchingParameters_Auto(SUNNonlinearSolver NLS, sunrealtype newt_to_fp_threshold, long int newt_to_fp_delay, sunrealtype fp_to_newt_threshold, long int fp_to_newt_delay)

This function sets the parameters that control the switching behavior of the SUNNonlinearSolver_Auto module.

Parameters:

NLS – the nonlinear solver object returned by SUNNonlinSol_Auto().
newt_to_fp_threshold – the threshold for switching from Newton to fixed-point (i.e., \(\beta\)). Nonnegative values are accepted; negative values reset the default 2.0.
newt_to_fp_delay – the minimum number of nonlinear solves after a switch before switching from Newton to fixed-point is allowed. Nonnegative values are accepted; negative values reset the default 10.
fp_to_newt_threshold – the threshold for switching from fixed-point to Newton (i.e., \(\alpha\)). Nonnegative values are accepted; negative values reset the default 0.8.
fp_to_newt_delay – the minimum number of nonlinear solves after a switch before switching from fixed-point to Newton is allowed. Nonnegative values are accepted; negative values reset the default 1.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

Note

If supported by the SUNNonlinearSolver implementation, this routine will be called by SUNNonlinSolSetOptions() when using the key NLSid.switching_parameters followed by four values in the order newt_to_fp_threshold, newt_to_fp_delay, fp_to_newt_threshold, and fp_to_newt_delay.

SUNErrCode SUNNonlinSolGetFixedPointSolver_Auto(SUNNonlinearSolver NLS, SUNNonlinearSolver *fp_nls)

This function returns the underlying fixed-point solver so that users may configure fixed-point-specific options directly, e.g., SUNNonlinSolSetMaxIters().

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
fp_nls – a pointer to the underlying fixed-point solver object.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

SUNErrCode SUNNonlinSolGetNewtonSolver_Auto(SUNNonlinearSolver NLS, SUNNonlinearSolver *newton_nls)

This function returns the underlying Newton solver so that users may configure Newton-specific options directly, e.g., SUNNonlinSolSetMaxIters().

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
newton_nls – a pointer to the underlying Newton solver object.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

SUNErrCode SUNNonlinSolGetActiveSolverType_Auto(SUNNonlinearSolver NLS, SUNNonlinSolAutoType *active_solver_type)

This function returns the currently active sub-solver.

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
active_solver_type – the currently active sub-solver type.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

SUNErrCode SUNNonlinSolGetSwitchCount_Auto(SUNNonlinearSolver NLS, long int *switch_count)

This function returns the number of automatic switches between the fixed-point and Newton sub-solvers over the life of the solver.

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
switch_count – the total number of automatic solver switches.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

SUNErrCode SUNNonlinSolGetTotalNumItersByType_Auto(SUNNonlinearSolver NLS, long int *fp_iters, long int *newt_iters)

This function returns the total nonlinear iteration counts accumulated by the fixed-point and Newton sub-solvers over the life of the solver.

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
fp_iters – the total number of fixed-point iterations.
newt_iters – the total number of Newton iterations.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

SUNErrCode SUNNonlinSolGetTotalNumConvFailsByType_Auto(SUNNonlinearSolver NLS, long int *fp_nconvfails, long int *newt_nconvfails)

This function returns the total nonlinear convergence failures accumulated by the fixed-point and Newton sub-solvers over the life of the solver.

Parameters:

NLS – a SUNNonlinearSolver object returned by SUNNonlinSol_Auto().
fp_nconvfails – the total number of fixed-point convergence failures.
newt_nconvfails – the total number of Newton convergence failures.

Returns:

SUN_SUCCESS if successful, otherwise an error code.

11.10. The SUNNonlinSol_PetscSNES implementation

This section describes the SUNNonlinSol interface to the PETSc SNES nonlinear solver(s). To enable the SUNonlinSol_PetscSNES module, SUNDIALS must be configured to use PETSc. Instructions on how to do this are given in §1.1.3.32. To access the SUNNonlinSol_PetscSNES module, include the header file sunnonlinsol/sunnonlinsol_petscsnes.h. The library to link to is libsundials_sunnonlinsolpetsc.lib where .lib is typically .so for shared libraries and .a for static libraries. Users of the SUNNonlinSol_PetscSNES module should also see §8.14 which discusses the NVECTOR interface to the PETSc Vec API.

11.10.1. SUNNonlinSol_PetscSNES description

The SUNNonlinSol_PetscSNES implementation allows users to utilize a PETSc SNES nonlinear solver to solve the nonlinear systems that arise in the SUNDIALS integrators. Since SNES uses the KSP linear solver interface underneath it, the SUNNonlinSol_PetscSNES implementation does not interface with SUNDIALS linear solvers. Instead, users should set nonlinear solver options, linear solver options, and preconditioner options through the PETSc SNES, KSP, and PC APIs.

Important usage notes for the SUNNonlinSol_PetscSNES implementation:

The SUNNonlinSol_PetscSNES implementation handles calling SNESSetFunction at construction. The actual residual function \(F(y)\) is set by the SUNDIALS integrator when the SUNNonlinSol_PetscSNES object is attached to it. Therefore, a user should not call SNESSetFunction on a SNES object that is being used with SUNNonlinSol_PetscSNES. For these reasons it is recommended, although not always necessary, that the user calls SUNNonlinSol_PetscSNES() with the new SNES object immediately after calling SNESCreate.
The number of nonlinear iterations is tracked by SUNDIALS separately from the count kept by SNES. As such, the function SUNNonlinSolGetNumIters() reports the cumulative number of iterations across the lifetime of the SUNNonlinearSolver object.
Some “converged” and “diverged” convergence reasons returned by SNES are treated as recoverable convergence failures by SUNDIALS. Therefore, the count of convergence failures returned by SUNNonlinSolGetNumConvFails() will reflect the number of recoverable convergence failures as determined by SUNDIALS, and may differ from the count returned by SNESGetNonlinearStepFailures.
The SUNNonlinSol_PetscSNES module is not currently compatible with the CVODES or IDAS staggered or simultaneous sensitivity strategies.

11.10.2. SUNNonlinearSolver_PetscSNES functions

The SUNNonlinSol_PetscSNES module provides the following constructor for creating a SUNNonlinearSolver object.

SUNNonlinearSolver SUNNonlinSol_PetscSNES(N_Vector y, SNES snes, SUNContext sunctx)

This creates a SUNNonlinSol object that wraps a PETSc SNES object for use with SUNDIALS. This will call SNESSetFunction on the provided SNES object.

Arguments:

snes – a PETSc SNES object.
y – a N_Vector object of type NVECTOR_PETSC that is used as a template for the residual vector.
sunctx – the SUNContext object (see §1.3)

Return value:

A SUNNonlinSol object if the constructor exits successfully, otherwise it will be NULL.

Warning

This function calls SNESSetFunction and will overwrite whatever function was previously set. Users should not call SNESSetFunction on the SNES object provided to the constructor.

The SUNNonlinSol_PetscSNES module implements all of the functions defined in §11.1.1–§11.1.3 except for SUNNonlinSolSetup(), SUNNonlinSolSetLSetupFn(), SUNNonlinSolSetLSolveFn(), SUNNonlinSolSetConvTestFn(), and SUNNonlinSolSetMaxIters().

The SUNNonlinSol_PetscSNES functions have the same names as those defined by the generic SUNNonlinSol API with _PetscSNES appended to the function name. Unless using the SUNNonlinSol_PetscSNES module as a standalone nonlinear solver the generic functions defined in §11.1.1–§11.1.3 should be called in favor of the SUNNonlinSol_PetscSNES specific implementations.

The SUNNonlinSol_PetscSNES module also defines the following user-callable functions.

SUNErrCode SUNNonlinSolGetSNES_PetscSNES(SUNNonlinearSolver NLS, SNES *snes)

This gets the SNES object that was wrapped.

Arguments:

NLS – a SUNNonlinSol object.
snes – a pointer to a PETSc SNES object that will be set upon return.

Return value:

A SUNErrCode

SUNErrCode SUNNonlinSolGetPetscError_PetscSNES(SUNNonlinearSolver NLS, PetscErrorCode *error)

This gets the last error code returned by the last internal call to a PETSc API function.

Arguments:

NLS – a SUNNonlinSol object.
error – a pointer to a PETSc error integer that will be set upon return.

Return value:

A SUNErrCode

SUNErrCode SUNNonlinSolGetSysFn_PetscSNES(SUNNonlinearSolver NLS, SUNNonlinSolSysFn *SysFn)

This returns the residual function that defines the nonlinear system.

Arguments:

NLS – a SUNNonlinSol object.
SysFn – the function defining the nonlinear system.

Return value:

A SUNErrCode

11.10.3. SUNNonlinearSolver_PetscSNES content

The content field of the SUNNonlinSol_PetscSNES module is the following structure.

struct _SUNNonlinearSolverContent_PetscSNES {
  int sysfn_last_err;
  PetscErrorCode petsc_last_err;
  long int nconvfails;
  long int nni;
  void *imem;
  SNES snes;
  Vec r;
  N_Vector y, f;
  SUNNonlinSolSysFn Sys;
};

These entries of the content field contain the following information:

sysfn_last_err – last error returned by the system defining function,
petsc_last_err – last error returned by PETSc,
nconvfails – number of nonlinear converge failures (recoverable or not),
nni – number of nonlinear iterations,
imem – SUNDIALS integrator memory,
snes – PETSc SNES object,
r – the nonlinear residual,
y – wrapper for PETSc vectors used in the system function,
f – wrapper for PETSc vectors used in the system function,
Sys – nonlinear system defining function.