"Action" functions in oomph-lib's black-box Newton solver

oomph-lib's Problem class provides a black-box Newton solver,
Problem::newton_solve(), which allows the driver code to be as compact as this:

main()
{
// Create the problem object
MyProblem problem;
// Solve the problem, using oomph-lib's default Newton solver
problem.newton_solve();
}

While the availability of a black-box Newton solver is helpful, it is sometimes necessary to interact with the solution process, e.g. to perform certain additional tasks during the Newton iteration, or to analyse the convergence (or divergence!) of the iteration in more detail. For this purpose oomph-lib's Problem class has a number of empty virtual member functions that are executed at key stages of the Newton iteration. Two of these functions (Problem::actions_before_newton_solve() and Problem::actions_after_newton_solve() ) are pure virtual functions that must be implemented in every specific Problem, though they may, of course, be implemented as empty functions. Other functions may be overloaded to provide more detailed control.


Problem::newton_solve()

The following table lists the main steps performed by oomph-lib's Newton solver and indicates at which point the various "action" functions are executed.

\[ \left. . \hspace{6cm} \right. \]

SEQUENCE OF MATHEMATICAL STEPS SEQUENCE OF FUNCTION CALLS / KEY CONTROL PARAMETERS
Step 1: Perform actions before solve.

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_before_newton_solve();
Step 2: Initialise counter for number of Newton iterations.

\[ c=0 \]

unsigned count=0;
Start loop over Newton iterations
Step 3: Perform actions before the next Newton step

\[ \left. . \hspace{6cm}\right. \]

Problem::actions_before_newton_step();
Only during the first Newton step:
Step 4a: Perform actions before Newton convergence check

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_before_newton_convergence_check();
Step 4b: Compute residual vector

\[ {\cal R}_i = {\cal R}_i\left(U_1,...,U_{N}\right) \mbox{\ \ \ \ for $i=1,...,N$} \]

Vector<double> residuals;
Problem::get_residuals(residuals);
Step 4c: Convergence check. If converged go to Step 13

\[ \mbox{If $\max_i |{\cal R}_i| < {\tt tol} \Longrightarrow $ converged.} \]

...
// Check convergence
if (max_res<Problem::Newton_solver_tolerance)
{
...
}
...
Step 5: Check if maximum residual or number of iterations exceed their maxima.

\[ \mbox{If \ \ } \left\{ \begin{array}{c} c > {\tt Max\_newton\_iterations} \\ \mbox{ or } \\ \max_i |{\cal R}_i| > {\tt Max\_residuals} \end{array} \right\} \Longrightarrow \mbox{diverged.} \]

...
// Check divergence/non-convergence
if ((max_res>Problem::Max_residuals) ||
(count == Problem::Max_newton_iterations))
{
...
// Throw an error
throw NewtonSolverError(count,max_res);
}
...
Step 6: Solve linear system for correction of unknowns

\[ \begin{array}{l} \mbox{Solve} \\ \\ \hspace{1cm}\sum_{j=1}^N {\cal J}_{ij} \ \delta U_j = {\cal R}_i \mbox{\ \ \ (for $i=1,...,N$)} \\ \\ \mbox{for $\delta U_j$ \ \ ($j=1,...,N$).} \\ \\ \mbox{Here the Jacobian matrix is given by} \\ \\ \hspace{1cm}{\cal J}_{ij} = \frac{\partial {\cal R}_i}{\partial U_j} \mbox{\ \ \ (for $i,j=1,...,N$)} \\ \end{array} \]

// Use the linear solver specified by
// Problem::linear_solver_pt() to assemble
// and solve the linear system.
...
Step 7: Apply the corrections to the solution

\[ U_i := U_i - \delta U_i \mbox{\ \ \ (for $i=1,...,N$)} \\ \]

// Update the unknowns
...
Step 8: Perform actions after Newton step

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_after_newton_step();
Step 9: Increment counter for number of Newton iterations

\[ c:=c+1 \]

count++;
Step 10: Perform actions before Newton convergence check

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_before_newton_convergence_check();
Step 11: Compute residual vector

\[ {\cal R}_i = {\cal R}_i\left(U_1,...,U_{N}\right) \mbox{\ \ \ \ for $i=1,...,N$} \]

Vector<double> residuals;
Problem::get_residuals(residuals);
Step 12: Convergence check. If converged go to Step 13, else go to Step 3.

\[ \mbox{If $\max_i |{\cal R}_i| < {\tt tol} \Longrightarrow$ converged.} \]

...
// Check convergence
if (max_res<Problem::Newton_solver_tolerance)
{
...
}
...
End loop over Newton iterations
Step 13: Converged. Perform actions after solve.

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_after_newton_solve();



The function Problem::actions_before_newton_solve() is often used to update boundary conditions when performing parameter studies in which these change; Problem::actions_after_newton_solve() may be used to automatically perform any post-processing when a solution has been obtained. The function Problem::actions_before_newton_convergence_check() makes it possible to update any "dependent" problem parameters, i.e. parameters that depend on one or more of the unknowns in the problem but are not updated automatically when the unknown is changed. This arises most frequently in free-boundary problems in which the position of the nodes in the "bulk" mesh is determined by an algebraic node update procedure. When the Newton method updates the Data values that determine the shape of the domain boundary in Step 7, the nodal positions themselves are not updated automatically unless the node update function is executed in Problem::actions_before_newton_convergence_check(). See the discussion of the free-boundary Poisson problem for an example of its use.


Optimisation for linear problems

Recall that, by default, oomph-lib regards all problems as nonlinear. Provided a good initial guess for the unknowns is available, the Newton solver will converge quadratically towards the solution. Within this framework, linear problems are simply special cases for which the Newton iteration converges in one step. However, if the problem is known to be linear, a few of the steps in the generic Newton iteration are unnecessary. For instance, it is not necessary to check the residual before or after the Newton step as we know that the exact solution will have been obtained (modulo roundoff errors) after step 7. The computation of the global residual vectors in steps 4 and 11 (which require a finite amount of cpu time) are therefore superfluous and are omitted if the "user" declares a problem as linear by setting the flag

bool Problem::Problem_is_nonlinear

which is initialised to true in the constructor of the Problem base class to false.

Problem::unsteady_newton_solve(...)

The following table lists the main steps performed by oomph-lib's unsteady Newton solver and indicates at which point the various "action" functions are executed. The arguments passed to the function are the value of the timestep dt and a boolean flag shift_values which indicates whether or not history values are to be shifted.

\[ \left. . \hspace{6cm} \right. \]

SEQUENCE OF MATHEMATICAL STEPS SEQUENCE OF FUNCTION CALLS / KEY CONTROL PARAMETERS
Step 1: Perform actions before implicit timestep.

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_before_implicit_timestep();
Step 2: Shift time values and dts, according to control flag.

\[ \left. . \hspace{6cm} \right. \]

if(shift_values) { shift_time_values(); }
Step 3: Advance global time and set current value of dt.

\[ t:=t+dt \]

time_pt()->time()+=dt;
time_pt()->dt()=dt;
Step 4: Loop over all timesteppers and set weights.

\[ \left. . \hspace{6cm} \right. \]

for(...)
{
time_stepper_pt(i)->set_weights();
}
Step 5: Solve the non-linear problem for this timestep with Newton's method.

\[ \left. . \hspace{6cm} \right. \]

Problem::newton_solve();
Step 6: Perform actions after implicit timestep.

\[ \left. . \hspace{6cm} \right. \]

Problem::actions_after_implicit_timestep();