# Numerical Solutions (e)

Like STOMP, eSTOMP uses Newton-Raphson linearization to solve the set of nonlinear algebraic equations formed by the discretized governing equations. This scheme has two major computational components: (1) to compute the Jacobian matrix and solution vector elements and (2) to solve the resulting linear system of equations.

## Newton-Raphson Iteration

The nonlinearities in the coupled flow and energy transport system of equations are resolved through the application of the iterative Newton-Raphson technique. The Newton-Raphson linearization technique is an iterative method for solving nonlinear algebraic equations of the form where f(x) is differentiable in x. The linearization concept approximates f(x) with suitable tangents.

#### NR Iteration (Conceptual Plot)

Each iteration yields a new estimate of x as the intersection of the tangent to the function f(x) at the previous estimate of x and the abscissa axis, according to In this formulation f(x) is considered the equation residual. For convergent systems, the residual decreases quadratically with iteration. In multiple variable systems, as with the coupled flow and energy transport system, the scalar function, f(x), is replaced with a vector function R(x), according to The vector function, R(x), represents the system of nonlinear algebraic equations produced from discretizing the conservation equations for component mass and energy. The vector of unknowns, x, represents the set of primary variables for the system, which are determined by the operational mode and phase conditions. This can be rewritten in terms of increments to the primary variables, according to The partial derivatives shown above form the Jacobian matrix.

Jacobian example for a 1D system

For simplification, a one-dimensional system involving the solution of the water mass, air mass, and energy conservation equations will be considered. The system of linear equations that result from applying the Newton-Raphson linearization technique to this system of nonlinear algebraic equations for a computational domain with “n” nodes is shown below. Here, each Jacobian matrix element represents a block matrix of order three, according to Each unknown vector element represents a vector of increments to the primary variables of order three, according to Each solution vector element represents a vector of equation residuals of order three, according to Info For a two-dimensional system involving the solution of the water mass, air mass, and energy conservation equations, the Jacobian matrix would contain two extra bands of block matrices: one below and one above the diagonal band. These extra bands would be located one half-band width from the main diagonal band, where the half-band width equaled the lesser of the number of nodes per row or column for a two-dimensional grid. A three-dimensional grid would contain four extra bands of block matrices, two below and two above the diagonal band. The furthest bands would be located one half-band width from the main diagonal band, where the halfband width equaled the least number of nodes in a plane. For example, a three-dimensional Cartesian grid with 20 nodes in the “X” coordinate direction, 30 nodes in the “Y” coordinate direction, and 40 nodes in the “Z” coordinate direction, would have a half-band width of 600.

# Linear System Solvers

The linear system of equations produced during Newton-Raphson iteration is solved using an iterative linear system solver. PETSc (The Portable, Extensible Toolkit for Scientific Computation) is the only option available for the linear solve for eSTOMP.  PETSc is a suite of data structures and routines that provide the building blocks for the implementation of large-scale application codes on parallel (and serial) computers.  Like SPLIB, PETSc offers several options for sparse linear system solvers. The user is referred to the PETSc manual or website for further information on each option. eSTOMP is currently compatible with version 3.1-p8. The PETSc library can be downloaded from the PETSc website.

PETSc uses the MPI standard for all message-passing communication which enables parallel execution on multiple processors. When eSTOMP is compiled, PETSc must already be compiled so that the eSTOMP executable can link with the solver library.

Defaults The default iterative solver is bi-conjugate gradient stabilized (BiCGSTAB/KSPBCGS) and the default precondition option is incomplete LU-zero factorization (ILU/PCILU). These options are recommended for most problems, but the user may view or change these options in the petsc.F source file. Alternatively, options can be set at run time (e.g., -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>), which will override those specified in the petsc.F source file.