Newton's method

Newton's method is intended to compute updates to some q by solving the following equation.

\begin{equation} \frac{\partial R(q)}{\partial q} \Delta q = -R(q) \end{equation}

In the most basic implementation of Newton's method in PDESolver, q corresponds to the solution, and f(q) corresponds to the residual evaluation of the currently selected physics module.

An example of a more sophisticated use of Newton's method is within the Crank-Nicolson timestepper, which adds another layer on top of the physics residual:

\begin{equation} \frac{\partial g(R(q))}{\partial q} \Delta q = -g(R(q)) \end{equation}

Features

PDESolver's Newton's method has a wide variety of features. It contains the Jacobian calculation routines, which can be performed currently using:

finite-differencing
complex-step

The Jacobian functions can act upon any arbitrary residual.

Additionally, the following matrix forms are supported:

Julia dense
Julia sparse
PETSc sparse
Matrix-free

The function that performs the Newton iteration is

NonlinearSolvers.newtonInner — Function.

This function contains the algorithm for Newton's method. It is intended to be used by other functions rather than invoked as a solver directly. For example newton uses newtonInner to solve steady problems while crank_nicolson uses it to solve the nonlinear problem arising from the time discretization.

For reference, Newton's method solves the problem R(q) = 0 using dR/dq delta q = -R(q) where R is the rsidual and q is the solution

On entry, eqn.q_vec must contain the initial guess for q. On exit, eqn.q_vec will contain the solution to f(q) = 0. eqn.q will also be consistent with eqn.q_vec, as will the send and receive buffers in eqn.shared_data

On exit, rhs_vec will have the residual corresponding to eqn.q_vec in it, with the imaginary part set to zero.

The AbstractLO and AbstractPC supplied inside the LinearSolver object must match the jac_type.

When doing inexact Newton-Krylov, newonInner modifies the tolerances of the linear solver. Users calling newtonInner repeatedly with the same linear solver object should reset the initial tolerances as needed.

Inputs:

newton_data: NewtonData object, typically obtained from setupNewton
mesh: a mesh object
sbp: an SBP operator
eqn: a solution data object
opts: options dictionary
ls: a LinearSolver with the preconditioner and linear operator fully initialized.
rhs_vec: vector to store R(q) in
ctx_residual: extra data required by rhs_func

The user must supply two functions, one to calculate the residual vector (referred to as rhs_vec), and another to compute the Jacobian.

rhs_func should compute (eqn.q_vec) -> (rhs_vec) and have the signature

rhs_func(mesh, sbp, eqn, opts, rhs_vec, ctx_residual, t=0.0)

Note that the ctx_residual passed into newtonInner is passed directly to rhs_func. The contents of ctx_residual can be whatever is needed by rhs_func to perform its computation. See physicsRhs for the specific requirements on rhs_func and for an example implementation.

The same ctx_residual passed into newtonInner is passed directly to calcLinearOperator..

This function supports jacobian/preconditioner freezing using the prefix "newton".

Aliasing restrictions: None. In particular, rhs_vec can alias eqn.res_vec, and this leads so some efficiency because it avoids needlessly copying data.

Newton's method

Features

Newton internals

NewtonData API

Linear Operators and Preconditioners