oomph-lib

"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.

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 "enslaved" 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.

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

Optimisation for linear problems.

PDF file